C     Last change:  ERB  10 Jul 2002    1:37 pm
       SUBROUTINE LAK3AL(ISUMRX,ISUMIR,LCCOND,ICLAKE,
     1 MXLKND,LKNODE,LCSTAG,IN,IOUT,ILKCB,NLAKES,INTRB,INDV,LCCNDF,
     2 LCLKPR,LCLKEV,ISTGLD,ISTGNW,IICS,IISUB,ISILL,LCWTDR,IFREFM,NROW,
     3 NCOL,NLAY,IBNLK,ILKBL,LKACC1,LKACC2,LKACC3,LKACC4,LKACC5,LKACC6,
     4 LKACC7,LKACC8,LKACC9,LKACC10,LKACC11,LKDRY,IBTMS,LKNCNT,LKKSUB,
     5 LKSADJ,LKFLXI,LKNCNS,LKSVT,LKJCLS,THETA,LCRNF,ITRSS,NSSITR,
     6 SSCNCR,LKSSMN,LKSSMX,LKNCN,LKDSR,LKCNN,LKCHN,IAREN,IUNITSFR,
     7 LSOVOL,NSS,IUNITGWT,LSLAKE,LSPPT,LSRNF,LSAUG,NSOL,IMSUB,IMSUB1,
     8 LSCGWL,LSSLAK,LSSWIN,LSSWOT,LSSPPT,LSCDRW,LSSRUN,LSGWIN,LSGWOT,
     9 LSLKSM,LSKLK,LSDONE,LSFLOB,LSRTCO,LSCLKO,LSALKI,LSALKO,ISTRIN,
     & ISTROT,LKLMRR,IDSTRT,
     * LKVI,ISTGLD2,LKCLKI,
     * LKCUM1,LKCUM2,LKCUM3,LKCUM4,LKCUM5,
     * LKCUM6,LKCUM7,LKCUM8,LKCUM9)
C
C----- USGS VERSION 3; AUG99 LAK3AL
C
C      ******************************************************************
C      ALLOCATE ARRAY STORAGE FOR LAKES
C      ******************************************************************
C
C      ------------------------------------------------------------------
C      SPECIFICATIONS:
C      ------------------------------------------------------------------
C
C1------IDENTIFY PACKAGE AND INITIALIZE LKNODE.
      WRITE(IOUT,1) IN
      LKNODE=0
C
C2------READ NLAKES, ILKCB.
C
      IF(IFREFM.EQ.0) THEN
         READ(IN,'(2I10)')NLAKES,ILKCB
         IF (ITRSS.EQ.0) THEN
            READ(IN,'(F10.2,I10,F10.2)') THETA,NSSITR,SSCNCR
         ELSE IF(ITRSS.GT.0) THEN
            READ(IN,'(F10.2)') THETA
            NSSITR=0
            SSCNCR=0.0
         ELSE
            WRITE(IOUT,14)
   14       FORMAT(1X,/,1X,
     1     'LAK CANNOT RUN WITH MIXED TRANSIENT AND STEADY STATE')
        call stopfile ! emrl
            STOP
         END IF
      ELSE
         READ(IN,*) NLAKES,ILKCB
         IF (ITRSS.EQ.0) THEN
            READ(IN,*) THETA,NSSITR,SSCNCR
         ELSE IF(ITRSS.GT.0) THEN
            READ(IN,*) THETA
            NSSITR=0
            SSCNCR=0.0
         ELSE
            WRITE(IOUT,14)
        call stopfile ! emrl
            STOP
         END IF
      END IF
C
C  VALUE OF MXLKND (NUMBER OF LAKE-AQUIFER INTERFACES) IS AN ESTIMATE.
C    TO SAVE MEMORY, REDUCE ITS SIZE IF APPROPRIATE.
C    IF MXLKND TOO SMALL, ERROR MESSAGE WILL BE PRINTED.
      MXLKND=NCOL*NROW*NLAY
C
      IF(NLAKES.LT.0) NLAKES=0
      WRITE(IOUT,5) MXLKND,NLAKES
      IF (ILKCB.GT.0) WRITE(IOUT,7) ILKCB
      IF (ILKCB.LE.0) WRITE(IOUT,9)
      WRITE(IOUT,22) THETA
      IF(ITRSS.EQ.0.AND.NSSITR.EQ.0) NSSITR = 50
      IF(ITRSS.EQ.0.AND.SSCNCR.EQ.0.0) SSCNCR = 0.01
      IF(ITRSS.EQ.0) WRITE(IOUT,23) NSSITR, SSCNCR
   23 FORMAT(/1X,'STEADY-STATE SOLUTION FOR LAKES.  MAXIMUM NUMBER OF IT
     1ERATIONS = ',I4,3X,'CONVERGENCE CRITERION = ',1PE9.2)
1     FORMAT(/1X,'LAK3 -- LAKE PACKAGE, VERSION 3, 6/28/99',
     1' INPUT READ FROM UNIT',I3)
5     FORMAT(1X,'SPACE ALLOCATION FOR',I7,' GRID CELL FACES ADJACENT TO
     1LAKES'/1X,'MAXIMUM NUMBER OF LAKES IS',I3, ' FOR THIS SIMULATION')
7     FORMAT(1X,'CELL-BY-CELL FLOWS WILL BE RECORDED ON UNIT',I5)
9     FORMAT(1X,'CELL-BY-CELL SEEPAGES WILL NOT BE PRINTED OR SAVED')
   22 FORMAT(/1X,'THETA = ',F10.2,'  METHOD FOR UPDATING LAKE STAGES IN
     1ITERATIONS OF THE SOLUTION FOR AQUIFER HEADS.'/20X,'0.0 IS EXPLICI
     2T, 0.5 IS CENTERED, AND 1.0 IS FULLY IMPLICIT.')
C3------SET LCCOND EQUAL TO ADDRESS OF FIRST UNUSED SPACE IN X.
200   LCCOND=ISUMRX
C
C4------CALCULATE AMOUNT OF SPACE NEEDED FOR CONDUCTIVITY
C       FACTOR AND LAKEBED CONDUCTANCE LIST.
C
      ISPA=MXLKND
      ISUMRX=ISUMRX+ISPA
      LCCNDF=ISUMRX
      ISPA1=MXLKND
      ISUMRX=ISUMRX+ISPA1
C
C5------CALCULATE AMOUNT OF SPACE NEEDED FOR ILAKE LIST.
      ICLAKE=ISUMIR
      ISPB=5*MXLKND
      ISUMIR=ISUMIR+ISPB
C
C6------CALCULATE AMOUNT OF SPACE NEEDED FOR NLAKES LIST (STAGES, PRCPLK,
C       EVAPLK, AND WTHDRW ARRAYS).
      LCSTAG=ISUMRX
      ISPC=NLAKES
      ISUMRX=ISUMRX+ISPC
      LCLKPR=ISUMRX
      ISPC1=NLAKES
      ISUMRX=ISUMRX+ISPC1
      LCLKEV=ISUMRX
      ISPC2=NLAKES
      ISUMRX=ISUMRX+ISPC2
      LCWTDR=ISUMRX
      ICR=NLAKES
      ISUMRX=ISUMRX+ICR
      LCRNF=ISUMRX
      IRNF=NLAKES
      ISUMRX=ISUMRX+IRNF
C
C7---  CALCULATE AMOUNT OF SPACE NEEDED FOR ITRB ARRAY.
      INTRB=ISUMIR
      IF(IUNITSFR.GT.0) THEN
        ISPD=NLAKES*NSS
        ISUMIR=ISUMIR+ISPD
      ELSE
        ISPD=0
      END IF
C
C8----  CALCULATE AMOUNT OF SPACE NEEDED FOR IDIV ARRAY.
      INDV=ISUMIR
      IF(IUNITSFR.GT.0) THEN
        ISPE=NLAKES*NSS
        ISUMIR=ISUMIR+ISPE
      ELSE
        ISPE=0
      END IF
C----  SET DUMMY VALUES FOR ISTRIN & ISTROT & IDSTRT IF IUNITSFR NOT ACTIVE
      IF (IUNITSFR.LE.0) THEN
         ISTRIN=ISUMRX
         ISTROT=ISUMRX
         IDSTRT=1
      END IF
C
C8----  CALCULATE AMOUNT OF SPACE NEEDED FOR STGOLD AND STGNEW.
      ISTGLD=ISUMRX
      ISPL=NLAKES
      ISUMRX=ISUMRX+ISPL
      ISTGNW=ISUMRX
      ISPM=NLAKES
      ISUMRX=ISUMRX+ISPM
C
C----  CALCULATE AMOUNT OF SPACE NEEDED FOR LAKE COALESCENSE PARAMETERS
C      ICS, ISUB, SILLVT, MSUB, MSUB1
      IICS=ISUMIR
      ICN=NLAKES
      ISUMIR=ISUMIR+ICN
      IISUB=ISUMIR
      ICO=NLAKES*NLAKES
      ISUMIR=ISUMIR+ICO
      ISILL=ISUMRX
      ICP=NLAKES*NLAKES
      ISUMRX=ISUMRX+ICP
      IMSUB=ISUMIR
      IMT=NLAKES*NLAKES
      ISUMIR=ISUMIR+IMT
      IMSUB1=ISUMIR
      IMT1=NLAKES
      ISUMIR=ISUMIR+IMT1
C
C---- CALCULATE AMOUNT OF SPACE NEEDED TO INPUT LAKE BORDER NODE AND BED
C     CONDUCTANCE INFORMATION
      IBNLK=ISUMIR
      ICS=NCOL*NROW*NLAY
      ISUMIR=ISUMIR + ICS
      ILKBL=ISUMRX
      ICT=NCOL*NROW*NLAY
      ISUMRX=ISUMRX+ICT
C
C---- PROVIDE SPACE FOR BUDGET ACCUMULATION TERMS
      LKACC1=ISUMRX
      IAC1=NLAKES
      ISUMRX=ISUMRX+IAC1
      LKACC2=ISUMRX
      IAC2=NLAKES
      ISUMRX=ISUMRX+IAC2
      LKACC3=ISUMRX
      IAC3=NLAKES
      ISUMRX=ISUMRX+IAC3
      LKACC4=ISUMRX
      IAC4=NLAKES
      ISUMRX=ISUMRX+IAC4
      LKACC5=ISUMRX
      IAC5=NLAKES
      ISUMRX=ISUMRX+IAC5
      LKACC6=ISUMRX
      IAC6=NLAKES
      ISUMRX=ISUMRX+IAC6
      LKACC7=ISUMRX
      IAC7=NLAKES
      ISUMRX=ISUMRX+IAC7
      LKACC8=ISUMRX
      IAC8=NLAKES
      ISUMRX=ISUMRX+IAC8
      LKACC9=ISUMRX
      IAC9=NLAKES
      ISUMRX=ISUMRX+IAC9
      LKACC10=ISUMRX
      IAC10=NLAKES
      ISUMRX=ISUMRX+IAC10
      LKACC11=ISUMRX
      IAC11=NLAKES
      ISUMRX=ISUMRX+IAC11
C
C---- PROVIDE SPACE FOR CUMULATIVE BUDGET TERMS
      LKCUM1=ISUMRX
      IAC1=NLAKES
      ISUMRX=ISUMRX+IAC1
      LKCUM2=ISUMRX
      IAC2=NLAKES
      ISUMRX=ISUMRX+IAC2
      LKCUM3=ISUMRX
      IAC3=NLAKES
      ISUMRX=ISUMRX+IAC3
      LKCUM4=ISUMRX
      IAC4=NLAKES
      ISUMRX=ISUMRX+IAC4
      LKCUM5=ISUMRX
      IAC5=NLAKES
      ISUMRX=ISUMRX+IAC5
      LKCUM6=ISUMRX
      IAC6=NLAKES
      ISUMRX=ISUMRX+IAC6
      LKCUM7=ISUMRX
      IAC7=NLAKES
      ISUMRX=ISUMRX+IAC7
      LKCUM8=ISUMRX
      IAC8=NLAKES
      ISUMRX=ISUMRX+IAC8
      LKCUM9=ISUMRX
      IAC9=NLAKES
      ISUMRX=ISUMRX+IAC9
C---  6 ARRAYS FOR LAK3BD
      LKNCNT=ISUMIR
      ISUMIR=ISUMIR+NLAKES
      LKKSUB=ISUMIR
      ISUMIR=ISUMIR+NLAKES
      LKSADJ=ISUMRX
      ISUMRX=ISUMRX+NLAKES
      LKFLXI=ISUMRX
      ISUMRX=ISUMRX+NLAKES
      LKNCNS=ISUMIR
      ISUMIR=ISUMIR+NLAKES
      LKSVT=ISUMRX
      ISUMRX=ISUMRX+NLAKES
      LKJCLS=ISUMIR
      ISUMIR=ISUMIR+NLAKES*NLAKES
      IAC12R=3*NLAKES
      IAC12I=3*NLAKES+NLAKES*NLAKES
C
C---- PROVIDE SPACE FOR LDRY ARRAY SHOWING WHICH LAKE CELLS HAVE GONE DRY
      LKDRY=ISUMIR
      ICDR=NROW*NCOL*NLAY
      ISUMIR=ISUMIR+ICDR
C
C---- PROVIDE SPACE FOR STORAGE OF LAKE BOTTOM ELEVATIONS
      IBTMS=ISUMRX
      IBT=NLAKES
      ISUMRX=ISUMRX+IBT
C
C---- PROVIDE SPACE FOR STORAGE OF MAXIMUM LAKE SURFACE AREAS
      IAREN=ISUMRX
      IAT=NLAKES
      ISUMRX=ISUMRX+IAT
C
C---- PROVIDE SPACE FOR STEADY-STATE CALCULATIONS
      LKNCN=ISUMIR
      ISS1=NLAKES
      ISUMIR=ISUMIR+ISS1
      LKDSR=ISUMRX
      ISS3=NLAKES
      ISUMRX=ISUMRX+ISS3
      LKSSMN=ISUMRX
      ISMN=NLAKES
      ISUMRX=ISUMRX+ISMN
      LKSSMX=ISUMRX
      ISMX=NLAKES
      ISUMRX=ISUMRX+ISMX
      LKLMRR=ISUMIR
      ISUMIR=ISUMIR+ISMX
C
C---- PROVIDE SPACE FOR IMPLICIT CALCULATIONS
      LKCNN=ISUMRX
      IIPT1=NLAKES
      ISUMRX=ISUMRX+IIPT1
      LKCHN=ISUMRX
      IIPT2=NLAKES
      ISUMRX=ISUMRX+IIPT2
C
C---- PROVIDE SPACE FOR GAGE COUNTERS
      LKVI=ISUMRX
      INVI=NLAKES
      ISUMRX=ISUMRX+INVI
      ISTGLD2=ISUMRX
      INSO=NLAKES
      ISUMRX=ISUMRX+INSO
      LKCLKI=ISUMRX
      ICLI=NLAKES*NSOL
      ISUMRX=ISUMRX+ICLI
C
C---- PROVIDE SPACE FOR LAKE-SOLUTE PARAMETERS
      LSLAKE=1
      LSAUG=1
      LSPPT=1
      LSRNF=1
      LSCGWL=1
      LSSLAK=1
      LSSWIN=1
      LSSWOT=1
      LSSPPT=1
      LSCDRW=1
      LSSRUN=1
      LSGWIN=1
      LSGWOT=1
      LSOVOL=1
      LSKLK=1
      LSDONE=1
      LSLKSM=1
      LSFLOB=1
      LSRTCO=1
      LSCLKO=1
      LSALKI=1
      LSALKO=1
      LSL1=0
      LSL2=0
      LSL3=0
      LSL4=0
      LSL5=0
      LSL6=0
      LSL7=0
      IF(IUNITGWT.GT.0) THEN
        LSLAKE=ISUMRX
        LSL1=NLAKES*NSOL
        ISUMRX=ISUMRX+LSL1
        LSAUG=ISUMRX
        LSL2=NLAKES*NSOL
        ISUMRX=ISUMRX+LSL2
        LSPPT=ISUMRX
        LSL3=NLAKES*NSOL
        ISUMRX=ISUMRX+LSL3
        LSRNF=ISUMRX
        ISUMRX=ISUMRX+LSL3
        LSCGWL=ISUMRX
        LSL5=NSOL
        ISUMRX=ISUMRX+LSL5
        LSSLAK=ISUMRX
        LSL6=NLAKES*NSOL
        ISUMRX=ISUMRX+LSL6
        LSSWIN=ISUMRX
        ISUMRX=ISUMRX+LSL6
        LSSWOT=ISUMRX
        ISUMRX=ISUMRX+LSL6
        LSSPPT=ISUMRX
        ISUMRX=ISUMRX+LSL6
        LSCDRW=ISUMRX
        ISUMRX=ISUMRX+LSL6
        LSSRUN=ISUMRX
        ISUMRX=ISUMRX+LSL6
        LSGWIN=ISUMRX
        ISUMRX=ISUMRX+LSL6
        LSGWOT=ISUMRX
        ISUMRX=ISUMRX+LSL6
        LSL7=NLAKES
        LSOVOL=ISUMRX
        ISUMRX=ISUMRX+LSL7
        LSLKSM=ISUMRX
        ISUMRX=ISUMRX+LSL5
        LSKLK=ISUMRX
        ISUMRX=ISUMRX+LSL7
        LSDONE=ISUMRX
        ISUMRX=ISUMRX+LSL7
        LSL4=MXLKND
        LSRTCO=ISUMRX
        ISUMRX=ISUMRX+LSL4
        LSFLOB=ISUMRX
        ISUMRX=ISUMRX+LSL4
        LSCLKO=ISUMRX
        ISUMRX=ISUMRX+LSL6
        LSALKI=ISUMRX
        ISUMRX=ISUMRX+LSL5
        LSALKO=ISUMRX
        ISUMRX=ISUMRX+LSL5
      END IF
C   ----------------------------------------------------------------
      ISPRX=ISPA+ISPC+ISPL
     1    +ISPM+ISPA1+ISPC1+ISPC2+ICP+ICR+ICT
     2    +IAC1+IAC2+IAC3+IAC4+IAC5+IAC6+IAC7+IAC8
     3    +IAC9+IAC10+IAC11+IAC12R+IBT+ISS3+IIPT1+IIPT2
     4    +IAT+LSL1+LSL2+2*LSL3+2*LSL4+4*LSL5+9*LSL6+3*LSL7
     5    +ISMN+2*ISMX+INVI+INSO+ICLI
     6    +IAC1+IAC2+IAC3+IAC4+IAC5+IAC6+IAC7+IAC8+IAC9
      ISPIR=ISPB+ISPD+ISPE+ICN+ICO+IMT+IMT1+ICS+ICDR+ISS1+IAC12I
C
C10------PRINT AMOUNT OF SPACED USED BY LAKE PACKAGE.
      WRITE (IOUT,18) ISPRX
18    FORMAT(1X,I10,' ELEMENTS IN RX ARRAY ARE USED BY LAK')
      WRITE (IOUT,19) ISPIR
19    FORMAT(1X,I10,' ELEMENTS IN IR ARRAY ARE USED BY LAK')
C
C11-----RETURN.
      IF(IUNITSFR.LE.0) NSS=1
      RETURN
      END
C
      SUBROUTINE LAK3RP(ILAKE,LKNODE,MXLKND,IN,IOUT,NLAKES,STAGES,
     1 PRCPLK,EVAPLK,BEDLAK,NTRB,NDV,ITRB,IDIV,KKPER,
     2 DELR,DELC,NCOL,NROW,NLAY,ICS,VOL,BOTM,NBOTM,
     3 ISUB,SILLVT,ICMX,NCLS,WTHDRW,LWRT,IFREFM,LKARR,BDLKNC,LKARR1,
     4 BDLKN1,NODES,BOTTMS,RNF,BGAREA,IUNITSFR,NSS,IUNITGWT,CLAKE,
     5 CAUG,CPPT,CRNF,NSOL,IOUTS,SSMN,SSMX,ISS,VOLINIT,CLAKINIT,
     6 CUMPPT,CUMEVP,CUMRNF,
     1 CUMGWI,CUMGWO,CUMSWI,CUMSWO,CUMWDR,CUMFLX)
C
C----- USGS VERSION 3;  AUG99 LAK3RP
C
C     ******************************************************************
C       READ INPUT DATA FOR THE LAKE PACKAGE.
C     ------------------------------------------------------------------
C     SPECIFICATIONS:
C     ------------------------------------------------------------------
      COMMON /DISCOM/LBOTM(200),LAYCBD(200)
      CHARACTER*24 ANAME(2)
      CHARACTER*30 LFRMAT
      DIMENSION BEDLAK(MXLKND),ILAKE(5,MXLKND)
      DIMENSION STAGES(NLAKES),PRCPLK(NLAKES),EVAPLK(NLAKES)
      DIMENSION DELR(NCOL),RNF(NLAKES),DELC(NROW)
      DIMENSION ITRB(NLAKES,NSS),IDIV(NLAKES,NSS),WTHDRW(NLAKES)
      DIMENSION ICS(NLAKES),ISUB(NLAKES,NLAKES),SILLVT(NLAKES,NLAKES),
     1  LKARR(NODES),BDLKNC(NODES),LKARR1(NCOL,NROW,NLAY),
     2  BDLKN1(NCOL,NROW,NLAY),BOTTMS(NLAKES),BOTM(NCOL,NROW,0:NBOTM),
     3  BGAREA(NLAKES),CLAKE(NLAKES,NSOL),CAUG(NLAKES,NSOL),
     4  CPPT(NLAKES,NSOL),CRNF(NLAKES,NSOL),VOL(NLAKES),
     5  SSMN(NLAKES),SSMX(NLAKES)
      DIMENSION VOLINIT(NLAKES),CLAKINIT(NLAKES,NSOL)
      DIMENSION CUMPPT(NLAKES),CUMEVP(NLAKES),CUMRNF(NLAKES),
     1 CUMGWI(NLAKES),CUMGWO(NLAKES),CUMSWI(NLAKES),CUMSWO(NLAKES),
     2 CUMWDR(NLAKES),CUMFLX(NLAKES)
      DATA ANAME(1)/'           LAKE ID ARRAY'/
      DATA ANAME(2)/'  LAKEBED LEAKANCE ARRAY'/
C
C     ------------------------------------------------------------------
C1A-----IF MXLKND IS LESS THAN 1, THEN LAKE IS INACTIVE. RETURN.
      IF(MXLKND.LT.1) RETURN
C
C1A1----READ INITIAL CONDITIONS FOR ALL LAKES (ONLY READ ONCE)
      IF (KKPER.EQ.1) THEN
         IF(ISS.NE.0) WRITE (IOUT,20)
         IF(ISS.EQ.0) WRITE (IOUT,820)
         IF (IUNITGWT.EQ.0) THEN
            DO 30 LM=1,NLAKES
               IF (IFREFM.EQ.0) THEN
                  IF(ISS.NE.0) READ (IN,'(3F10.4)') STAGES(LM),SSMN(LM),
     1              SSMX(LM)
                  IF(ISS.EQ.0) READ (IN,'(3F10.4)') STAGES(LM)
               ELSE
                  IF(ISS.NE.0) READ (IN,*) STAGES(LM),SSMN(LM),SSMX(LM)
                  IF(ISS.EQ.0) READ (IN,*) STAGES(LM)
               END IF
            IF(ISS.NE.0) WRITE (IOUT,22) LM,STAGES(LM),SSMN(LM),SSMX(LM)
            IF(ISS.EQ.0) WRITE (IOUT,22) LM,STAGES(LM)
 30         CONTINUE
         ELSE
            WRITE (IOUTS,21) NSOL
            WRITE (LFRMAT,23) NSOL
            DO 35 LM=1,NLAKES
               IF (IFREFM.EQ.0) THEN
                 IF(ISS.NE.0) READ(IN,'(100F10.4)') STAGES(LM),SSMN(LM),
     1              SSMX(LM),(CLAKE(LM,ISOL),ISOL=1,NSOL)
                 IF(ISS.EQ.0) READ (IN,'(100F10.4)') STAGES(LM),
     1                        (CLAKE(LM,ISOL),ISOL=1,NSOL)
               ELSE
                 IF(ISS.NE.0) READ (IN,*) STAGES(LM),SSMN(LM),SSMX(LM),
     1                        (CLAKE(LM,ISOL),ISOL=1,NSOL)
                 IF(ISS.EQ.0) READ (IN,*) STAGES(LM),
     1                        (CLAKE(LM,ISOL),ISOL=1,NSOL)
               END IF
            IF(ISS.NE.0) WRITE (IOUT,22) LM,STAGES(LM),SSMN(LM),SSMX(LM)
            IF(ISS.EQ.0) WRITE (IOUT,22) LM,STAGES(LM)
 35           WRITE (IOUTS,LFRMAT) LM,(CLAKE(LM,ISOL),ISOL=1,NSOL)
cgage
            CLAKINIT=CLAKE
         END IF
      END IF
      WRITE (IOUT,'(/)')
 20   FORMAT(///1X,'INITIAL LAKE STAGE:  LAKE    STAGE    SS MIN    SS M
     1AX'/)
 820  FORMAT (/1X,'INITIAL LAKE STAGE:  LAKE    STAGE'/)
 21   FORMAT (//1X,'INITIAL LAKE CONCENTRATIONS:  LAKE   CONCENTRATION (
     1NSOL =',I3,')'/)
 22   FORMAT (22X,I3,3F10.3)
 23   FORMAT ('(31X,I3,3X,1P',I3,'(E12.3))')
C
C1B-----READ ITMP (FLAG TO REUSE LAKE-GEOMETRY DATA).
      IF(IFREFM.EQ.0) THEN
         READ(IN,'(3I10)') ITMP, ITMP1, LWRT
      ELSE
         READ(IN,*) ITMP, ITMP1, LWRT
      END IF
C
C2A-----IF ITMP < 0 THEN REUSE LAKE CONFIGURATION DATA FROM LAST STRESS
C       PERIOD.
      IF(ITMP.GE.0) GO TO 50
      WRITE (IOUT,'(/)')
      WRITE(IOUT,2)
    2 FORMAT(1H ,'REUSING LAKE CONFIGURATION DATA FROM LAST STRESS PERIO
     1D'/)
      GO TO 800
C
C4------IF THERE ARE NO LAKE NODES THEN RETURN.
   50 LKNODE = 0
      IF(ITMP.EQ.0) RETURN
C
C   INITIALIZE BGAREA
      DO 60 LK=1,NLAKES
      BGAREA(LK)=0.0
   60 CONTINUE
C
C5------READ INTEGER ARRAYS THAT DEFINE THE POSITIONS OF ALL LAKES IN
C5A     EACH MODEL GRID LAYER.  THEN READ ARRAYS OF LAKEBED CONDUCTANCES
C5B     IN EACH LAYER.
C
C   READ ARRAY OF LAKE ID'S, LAYER BY LAYER
C
      NCR = NCOL*NROW
      DO 125 K=1,NLAY
      KK = K
      LOC = 1 + (K-1)*NCR
      CALL U2DINT(LKARR(LOC),ANAME(1),NROW,NCOL,KK,IN,IOUT)
  125 CONTINUE
C
C   CHECK THAT ALL ENTRIES ARE VALID LAKE ID NUMBERS OR ZERO
C
      DO 130 K=1,NLAY
      DO 130 I=1,NCOL
      DO 130 J=1,NROW
      IF(LKARR1(I,J,K).GT.0.AND.LKARR1(I,J,K).LE.NLAKES) GO TO 130
      LKARR1(I,J,K)=0
  130 CONTINUE
C
C   READ ARRAY OF BED LEAKANCES, LAYER BY LAYER
C
      WRITE (IOUT,'(/)')
      DO 135 K=1,NLAY
      KK = K
      LOC = 1 + (K-1)*NCR
      CALL U2DREL(BDLKNC(LOC),ANAME(2),NROW,NCOL,KK,IN,IOUT)
  135 CONTINUE
C
        WRITE(IOUT,36)
        WRITE(IOUT,4)
36    FORMAT(/7X,'LOCATIONS, LAKE #, INTERFACE TYPE FOR GRID CELLS',
     1 ' ADJACENT TO LAKES:',5X,/
     3 5X,71('-'))
4     FORMAT(5X,'LAYER #',4X,'ROW #',4X,'COLUMN #',3X,'LAKE #',
     1       2X,'INTERFACE TYPE',2X,'LAKEBED LEAKANCE')
C
C   IDENTIFY LAKE BORDER CELLS, ASSIGN CELL TYPE ID'S, COMPUTE AND
C     ASSIGN LAKE-AQUIFER INTERFACE CONDUCTANCES.
C
      M = 0
      DO 180 I=1,NCOL
      DO 180 J=1,NROW
      K = 1
      IF(LKARR1(I,J,K).EQ.0) GO TO 150
      IF(NLAY.EQ.1) GO TO 145
C   Keep searching in vertical direction until non-lake cell is found, and define
C     interface there ("K" for interface is layer below bottom of lake)
      DO 140 K=2,NLAY
      IF(LKARR1(I,J,K).EQ.0) GO TO 145
  140 CONTINUE
      GO TO 145
C
C   VERTICAL LAKEBED INTERFACE (TYPE 0) DETECTED
C
  145 M = M + 1
      IF(M.LE.MXLKND) GO TO 147
      WRITE(IOUT,149) I,J,K
  149 FORMAT(/1X,'MAXIMUM NUMBER OF GRID CELLS ADJACENT TO LAKES HAS BEE
     1N EXCEEDED WITH CELL ',3I5,'  REDEFINE VARIABLE MXLKND TO A LARGER
     2 VALUE IN MODULE LAK3AL')
        call stopfile ! emrl
      STOP
  147 ILAKE(1,M) = K
      ILAKE(2,M) = J
      ILAKE(3,M) = I
      IF(K.GT.1.AND.LKARR1(I,J,K).EQ.0) LID = LKARR1(I,J,K-1)
      IF(LKARR1(I,J,K).NE.0) LID = LKARR1(I,J,K)
      ILAKE(4,M) = LID
      ILAKE(5,M) = 6
      BEDLAK(M) = BDLKN1(I,J,K-1)
      IF(K.EQ.NLAY.AND.LKARR1(I,J,K).NE.0) BEDLAK(M) = 0.0
      BGAREA(LID) = BGAREA(LID) + DELC(J)*DELR(I)
      WRITE(IOUT,5) (ILAKE(I1,M),I1=1,5), BEDLAK(M)
5     FORMAT(5I10,10X,F10.5)
      IF(LKARR1(I,J,K).NE.0) GO TO 180
C
C   SEARCH FOR CELL(S) ADJACENT TO LAKE
C
  150 K2 = K
      DO 175 K1=K2,NLAY
      IF(I.NE.1) GO TO 1151
      IF(LKARR1(I+1,J,K1).EQ.0) GO TO 165
      GO TO 1153
 1151 IF(I.NE.NCOL) GO TO 1152
      IF(LKARR1(I-1,J,K1).EQ.0) GO TO 165
      GO TO 1153
 1152 IF(LKARR1(I+1,J,K1).EQ.0.AND.LKARR1(I-1,J,K1).EQ.0) GO TO 165
C
C   CELL(S) LATERALLY ADJACENT TO LAKE IN X-DIRECTION (TYPE 1) DETECTED
C
 1153 DO 160 N=1,2
      IF(N.EQ.2) GO TO 155
      IF(I.EQ.1) GO TO 160
      IF(LKARR1(I-1,J,K1).EQ.0) GO TO 160
      I2 = I-1
      IFACE=1
      GO TO 157
  155 IF(I.EQ.NCOL) GO TO 160
      IF(LKARR1(I+1,J,K1).EQ.0) GO TO 160
      I2 = I + 1
      IFACE=2
  157 M = M + 1
      IF(M.LE.MXLKND) GO TO 158
      WRITE(IOUT,149) I,J,K1
        call stopfile ! emrl
      STOP
  158 ILAKE(1,M) = K1
      ILAKE(2,M) = J
      ILAKE(3,M) = I
      ILAKE(4,M) = LKARR1(I2,J,K1)
      ILAKE(5,M) = IFACE
      BEDLAK(M) = BDLKN1(I,J,K1)
      K4 = K1 - 1
      DO 3158 K3=1,K4
      IF(LKARR1(I,J,K3).EQ.0) GO TO 3158
      GO TO 3162
 3158 CONTINUE
      BEDLAK(M) = BDLKN1(I,J,1)
 3162 CONTINUE
      WRITE(IOUT,5) (ILAKE(I1,M),I1=1,5), BEDLAK(M)
  160 CONTINUE
  165 IF(J.NE.1) GO TO 1161
      IF(LKARR1(I,J+1,K1).EQ.0) GO TO 175
      GO TO 1163
 1161 IF(J.NE.NROW) GO TO 1162
      IF(LKARR1(I,J-1,K1).EQ.0) GO TO 175
      GO TO 1163
 1162 IF(LKARR1(I,J+1,K1).EQ.0.AND.LKARR1(I,J-1,K1).EQ.0) GO TO 175
C
C   CELL(S) LATERALLY ADJACENT TO LAKE IN Y-DIRECTION (TYPE 2) DETECTED
C
 1163 DO 170 N=1,2
      IF(N.EQ.2) GO TO 172
      IF(J.EQ.1) GO TO 170
      IF(LKARR1(I,J-1,K1).EQ.0) GO TO 170
      J2 = J - 1
      IFACE=4
      GO TO 174
  172 IF(J.EQ.NROW) GO TO 170
      IF(LKARR1(I,J+1,K1).EQ.0) GO TO 170
      J2 = J + 1
      IFACE=3
  174 M = M + 1
      IF(M.LE.MXLKND) GO TO 176
      WRITE(IOUT,149) I,J,K1
        call stopfile ! emrl
      STOP
  176 ILAKE(1,M) = K1
      ILAKE(2,M) = J
      ILAKE(3,M) = I
      ILAKE(4,M) = LKARR1(I,J2,K1)
      ILAKE(5,M) = IFACE
      BEDLAK(M) = BDLKN1(I,J,K1)
      K4 = K1 - 1
      DO 4158 K3=1,K4
      IF(LKARR1(I,J,K3).EQ.0) GO TO 4158
      GO TO 4162
 4158 CONTINUE
      BEDLAK(M) = BDLKN1(I,J,1)
 4162 CONTINUE
      WRITE(IOUT,5) (ILAKE(I1,M),I1=1,5), BEDLAK(M)
  170 CONTINUE
  175 CONTINUE
  180 CONTINUE
      WRITE(IOUT,195) M
  195 FORMAT(/5X,'NUMBER OF LAKE-AQUIFER CELL INTERFACES = ',I5)
      LKNODE = M
C
C   SET LAKE BOTTOM ELEVATIONS
      DO 295 LK=1,NLAKES
  295 BOTTMS(LK) = 999999
C
      DO 350 II=1,LKNODE
      K = ILAKE(1,II)
      J = ILAKE(2,II)
      I = ILAKE(3,II)
C  Convert ILAKE(5,II):  1 and 2 are type 1,  3 and 4 are type 2, 6 is type 0
      NTYP = (ILAKE(5,II)+1)/2
      IF(NTYP.EQ.3) NTYP=0
      IF(NTYP.EQ.0) THEN
        LAKE = ILAKE(4,II)
        IF(K.GT.1) BOTLK = BOTM(I,J,LBOTM(K-1))
        IF(K.EQ.NLAY.AND.LKARR1(I,J,K).GT.0) BOTLK = BOTM(I,J,LBOTM(K))
        IF(BOTLK.LT.BOTTMS(LAKE)) BOTTMS(LAKE) = BOTLK
      END IF
  350 CONTINUE
C
C-- COMPUTE AND PRINT STAGE/VOLUME TABLES WHEN MORE THAN ONE LAYER
C
      IF(NLAY.EQ.1) GO TO 1331
      DO 1330 L1=1,NLAKES
      WRITE(IOUT,1306) L1
 1306 FORMAT(//1X,'STAGE/VOLUME RELATION FOR LAKE',I3//3X,'STAGE',
     1        4X,'VOLUME'/)
      EVOL = 0.0
      GTSDPH = 0.0
      TOPMST = BOTTMS(L1)
      TBELV = BOTTMS(L1)
      DO 1340 I=1,NCOL
      DO 1340 J=1,NROW
      IF(LKARR1(I,J,1).NE.L1) GO TO 1340
      IF(BOTM(I,J,LBOTM(1)).GT.TOPMST) TOPMST = BOTM(I,J,LBOTM(1))
      DTHK = BOTM(I,J,LBOTM(1)) - BOTM(I,J,LBOTM(2))
      IF(DTHK.GT.GTSDPH) GTSDPH = DTHK
 1340 CONTINUE
      TOPMST = TOPMST + GTSDPH
      TBNC = (TOPMST-BOTTMS(L1))/100.0
      WRITE(IOUT,1315) TBELV, EVOL
 1315 FORMAT(F10.4,1X,1PE10.3)
      DO 1325 K=1,150
      EVOL=0.0
      TBELV = TBELV + TBNC
      DO 1320 I=1,NCOL
      DO 1320 J=1,NROW
      IF(LKARR1(I,J,1).NE.L1) GO TO 1320
      DO 1318 K2=1,NLAY
      IF(LKARR1(I,J,K2).EQ.0) GO TO 1319
 1318 CONTINUE
      BOTIJ = BOTM(I,J,LBOTM(NLAY))
      GO TO 1313
 1319 BOTIJ = BOTM(I,J,LBOTM(K2-1))
 1313 IF(TBELV.LE.BOTIJ) GO TO 1320
      DV = (TBELV-BOTIJ)*DELC(J)*DELR(I)
      EVOL = EVOL + DV
 1320 CONTINUE
      WRITE(IOUT,1315) TBELV, EVOL
 1325 CONTINUE
      WRITE(IOUT,1326)
 1326 FORMAT(120X)
 1330 CONTINUE
 1331 CONTINUE
C
C-- INITIALIZE STREAM INFLOW SEGMENT ARRAY TO ZERO.
      DO 400 LNUM=1,NLAKES
         DO 400 LNP=1,NSS
400         ITRB(LNUM,LNP)=0
C
C-- INITIALIZE STREAM OUTFLOW SEGMENT ARRAY TO ZERO.
      DO 500 LNUM=1,NLAKES
         DO 500 LNP=1,NSS
500         IDIV(LNUM,LNP)=0
C
      IF(IUNITSFR.LE.0) THEN
         NDV=0
         NTRB=0
      END IF
C
C
C--  READ LINKAGE PARAMETERS FOR COALESCING LAKES
C
C    FOR EACH CONNECTED LAKE SYSTEM, READ LAKE NUMBERS OF CENTER LAKES
C    AND ADJOINING LAKES AND SILL ELEVATIONS.  ENTER CARD IMAGES
C    FOR SUBLAKE SYSTEMS EVEN IF LINKED TO MAIN LAKE SYSTEM.  SYSTEMS
C    MUST BE ORDERED HEIRARCHICALLY.
C
      ICMX = 0
      NCLS=0
      IF(IFREFM.EQ.0) THEN
        READ(IN,'(I5)') NSLMS
      ELSE
        READ(IN,*) NSLMS
      END IF
      WRITE(IOUT,680) NSLMS
  680 FORMAT(/1X,'NUMBER OF CONNECTED LAKE SYSTEMS IN SIMULATION IS ',I3
     1)
      IF(NSLMS.LE.0) GO TO 760
      DO 700 IS=1,NSLMS
      IF(IFREFM.EQ.0) THEN
        READ(IN,'(16I5)',END=750) IC,(ISUB(IS,I),I=1,IC)
      ELSE
        READ(IN,*,END=750) IC,(ISUB(IS,I),I=1,IC)
      END IF
      IF(IC.LE.0) GO TO 750
      IF(IC.GT.ICMX) ICMX=IC
      ICS(IS)=IC
      IC1 = IC - 1
      IF(IFREFM.EQ.0) THEN
        READ(IN,'(100F10.2)') (SILLVT(IS,I),I=1,IC1)
      ELSE
        READ(IN,*) (SILLVT(IS,I),I=1,IC1)
      END IF
      WRITE(IOUT,18) IS, ICS(IS), ISUB(IS,1)
   18 FORMAT(/10X,'SYSTEM',I3//2X,'NUMBER OF LAKES IN SYSTEM',I5,
     1  '  CENTER LAKE NUMBER',I5//1X,'SUBLAKE NUMBER',3X,
     2  'SILL ELEVATION'/)
      DO 715 JK=2,IC
  715 WRITE(IOUT,717) ISUB(IS,JK), SILLVT(IS,JK-1)
  717 FORMAT(8X,I2,8X,F10.2)
  700 CONTINUE
  750 CONTINUE
      NCLS=IS-1
      WRITE(IOUT,751) NCLS
  751 FORMAT(/1X,'READ DATA FOR',I5,' LAKE SYSTEMS'/)
  760 CONTINUE
C
C----- READ LAKE PRECIPITATION, EVAPORATION, RUNOFF, AND WITHDRAWAL RATES.
C      IF ITMP1 LT 0, SPECIFICATIONS FROM LAST STRESS PERIOD ARE USED.
C
  800 IF(ITMP1.GE.0) GO TO 801
      WRITE(IOUT,802)
  802 FORMAT(1H0,'REUSING RECH,ET,WITHDRAWAL RATES FROM LAST STRESS PERI
     1OD'/)
      RETURN
  801 WRITE(IOUT,7)
7     FORMAT(/1X,'LAKE',7X,'PRECIP',5X,'EVAP',5X,'RUNOFF',
     2     3X,'WITHDRAW',3X,'BOTTOM',5X,'AREA'/70('-'))
      IF (IUNITGWT.GT.0) WRITE (IOUTS,8)
 8    FORMAT (//1X,'LAKE',4X,'SOLUTE',6X,'CPPT',6X,'CRNF',6X,'CAUG'/)
      DO 300 LM=1,NLAKES
      IF(IFREFM.EQ.0) THEN
        READ(IN,'(4F10.4)') PRCPLK(LM),EVAPLK(LM),RNF(LM),WTHDRW(LM)
      ELSE
        READ(IN,*) PRCPLK(LM),EVAPLK(LM),RNF(LM),WTHDRW(LM)
      END IF
      WRITE(IOUT,9) LM,PRCPLK(LM),EVAPLK(LM),RNF(LM),WTHDRW(LM),
     1 BOTTMS(LM),BGAREA(LM)
9     FORMAT(1X,I3,4X,1P,3E10.3,1X,3E10.3)
      IF(IUNITGWT.LE.0) GO TO 300
      DO 850 ISOL=1,NSOL
      IF(IFREFM.EQ.0) THEN
         IF(WTHDRW(LM).LT.0.0) THEN
            READ(IN,'(3F10.4)')CPPT(LM,ISOL),CRNF(LM,ISOL),CAUG(LM,ISOL)
         ELSE
            READ(IN,'(2F10.4)')CPPT(LM,ISOL),CRNF(LM,ISOL)
         END IF
      ELSE
         IF(WTHDRW(LM).LT.0.0) THEN
            READ(IN,*) CPPT(LM,ISOL),CRNF(LM,ISOL),CAUG(LM,ISOL)
         ELSE
            READ(IN,*) CPPT(LM,ISOL),CRNF(LM,ISOL)
         END IF
      END IF
      IF(WTHDRW(LM).LT.0.0)
     1WRITE(IOUTS,840) LM,ISOL,CPPT(LM,ISOL),CRNF(LM,ISOL),CAUG(LM,ISOL)
      IF(WTHDRW(LM).GE.0.0)
     1WRITE(IOUTS,841) LM,ISOL,CPPT(LM,ISOL),CRNF(LM,ISOL)
  840 FORMAT(1X,I3,6X,I3,4X,1P,3E10.2)
  841 FORMAT(1X,I3,6X,I3,4X,1P,2E10.2)
  850 CONTINUE
      WRITE (IOUTS,'(/)')
  300 CONTINUE
      WRITE (IOUT,'(/)')
C
C--  Define Initial Lake Volume & Initialize Cumulative Budget Terms
      IF(KKPER.EQ.1) THEN
         DO 8400 LK=1,NLAKES
 8400    VOL(LK)=0.0
         DO 8450 LK=1,NLAKES
             CUMPPT(NLAKES)=0.0
             CUMEVP(NLAKES)=0.0
             CUMRNF(NLAKES)=0.0
             CUMGWI(NLAKES)=0.0
             CUMGWO(NLAKES)=0.0
             CUMSWI(NLAKES)=0.0
             CUMSWO(NLAKES)=0.0
             CUMWDR(NLAKES)=0.0
             CUMFLX(NLAKES)=0.0
 8450    CONTINUE
            DO 8900 L=1,LKNODE
               IL=ILAKE(1,L)
               IR=ILAKE(2,L)
               IC=ILAKE(3,L)
               LAKE=ILAKE(4,L)
C  Convert ILAKE(5,L):  1 and 2 are type 1,  3 and 4 are type 2, 6 is type 0
               ITYPE = (ILAKE(5,L)+1)/2
               IF(ITYPE.EQ.3) ITYPE=0
               IF(ITYPE.NE.0) GO TO 8900
               IF(IL.GT.1) BOTLK = BOTM(IC,IR,LBOTM(IL-1))
               IF(IL.EQ.NLAY.AND.LKARR1(IC,IR,IL).GT.0)
     1            BOTLK = BOTM(IC,IR,LBOTM(IL))
               IF(STAGES(LAKE).GT.BOTLK) THEN
                  AREA = DELR(IC)*DELC(IR)
                  VOL(LAKE)=VOL(LAKE)+(STAGES(LAKE)-BOTLK)*AREA
cgage
                  VOLINIT(LAKE)=VOL(LAKE)
               ENDIF
 8900       CONTINUE
      ENDIF
C
C7------RETURN
      RETURN
      END
C
      SUBROUTINE LAK3AD(KKPER,KKSTP,NLAKES,STGOLD,STGNEW,STAGES,
     1        NROW,NCOL,NLAY,FLOB,IUNITGWT,LKNODE,STGOLD2)
C
C-----VERSION 3 AUGUST 1999 LAK3AD
C
C     ******************************************************************
C     ADVANCE TO NEXT TIME STEP FOR TRANSIENT LAKE SIMULATION, AND COPY
C             INITIAL LAKE STAGES TO STGOLD FOR STEADY STATE.
C     ******************************************************************
C     SPECIFICATIONS:
C     ------------------------------------------------------------------
      DIMENSION STAGES(NLAKES),STGOLD(NLAKES),STGNEW(NLAKES)
      DIMENSION FLOB(LKNODE)
      DIMENSION STGOLD2(NLAKES)
C     ------------------------------------------------------------------
C
C1 --- COPY INITIAL LAKE STAGES TO STGOLD.
          DO 10 I=1,NLAKES
      IF(KKPER.EQ.1.AND.KKSTP.EQ.1) STGOLD(I)=STAGES(I)
      IF(KKPER.EQ.1.AND.KKSTP.EQ.1) STGOLD2(I)=STAGES(I)
10    IF(KKPER.EQ.1.AND.KKSTP.EQ.1) STGNEW(I)=STAGES(I)
C2 ----- IF NOT FIRST TIME STEP, OR FIRST STRESS PERIOD, UPDATE
C           STGOLD BY STGNEW.
      IF (KKPER.NE.1.OR.KKSTP.NE.1) THEN
            DO 30 K=1,NLAKES
               STGOLD(K)=STGNEW(K)
30             STGOLD2(K)=STGNEW(K)
      ENDIF
C
C-----Initialize FLOB array (stores cell by cell flux between lake and
C                            aquifer)
      IF (IUNITGWT.GT.0) THEN
        DO 50 LK=1,LKNODE
 50        FLOB(LK)=0.0
      END IF
C
C3------RETURN
      RETURN
      END
C
      SUBROUTINE LAK3ST(NFLG,NCOL,NROW,NLAY,IBOUND,LKNODE,ILAKE,MXLKND,
     1 NLAKES,STGOLD,BOTM,NBOTM)
C   ********************************************************************
C   SET IBOUND VALUES SO THAT RECHARGE AND EVAPOTRANSPIRATION (ET) WILL
C   BE ASSIGNED CORRECTLY UNDERNEATH DRYING LAKES (NFLG = 0), OR RESET
C   IBOUND AFTER RECHARGE AND ET ARE COMPUTED (NFLG = 1).
C   ********************************************************************
C
C   SPECIFICATIONS:
C
C-----------------------------------------------------------------------
      COMMON /DISCOM/LBOTM(200),LAYCBD(200)
      DIMENSION IBOUND(NCOL,NROW,NLAY),ILAKE(5,MXLKND),STGOLD(NLAKES),
     1 BOTM(NCOL,NROW,0:NBOTM)
C-----------------------------------------------------------------------
C
      IF(LKNODE.EQ.0) RETURN
      DO 10 L=1,LKNODE
C  Convert ILAKE(5,L):  1 and 2 are type 1,  3 and 4 are type 2, 6 is type 0
      ITYPE = (ILAKE(5,L)+1)/2
      IF(ITYPE.EQ.3) ITYPE=0
C
C-------ONLY CHANGE IBOUND FOR VERTICALLY ADJACENT NODE FACES
      IF(ITYPE.NE.0) GO TO 10
      IL = ILAKE(1,L)
      IR = ILAKE(2,L)
      IC = ILAKE(3,L)
C
C-------RESET AFTER EXECUTING RECHARGE OR ET ROUTINES
      IF(NFLG.EQ.1) GO TO 8
C
C-------RESET BEFORE EXECUTING RECHARGE OR ET ROUTINES
      IBOUND(IC,IR,IL-1) = -7
C
C-------THIS IS THE CORRECT ASSIGNMENT IF PORTION OF LAKE IN COLUMN
C       IS WET.
      LAKE = ILAKE(4,L)
      IF(STGOLD(LAKE).GT.BOTM(IC,IR,LBOTM(IL)-1)) GO TO 10
C
C-------IF PORTION OF LAKE IN NODE IS DRY, LET RECHARGE AND ET BE
C       APPLIED TO THE AQUIFER NODE UNDERNEATH THE LAKE BY SETTING
C       IBOUND EQUAL TO 0.
    8 IBOUND(IC,IR,IL-1) = 0
   10 CONTINUE
C
C3------RETURN
      RETURN
      END
C
      SUBROUTINE LAK3FM(LKNODE,MXLKND,ILAKE,HNEW,HCOF,RHS,
     1 IBOUND,NCOL,NROW,NLAY,NLAKES,STGOLD,CNDFCT,BOTM,NBOTM,IOUT,
     2 DELT,NSS,NTRB,NDV,ITRB,IDIV,STRIN,STROUT,STGNEW,WTHDRW,
     3 PRCPLK,EVAPLK,DELR,DELC,EVAP,PRECIP,SEEP,SURFA,SURFIN,SURFOT,
     4 THETA,RNF,KSTP,KITER,ISS,NSSITR,SSCNCR,SSMN,SSMX,DSTROT,
     5 NCNCVR,DSRFOT,SUMCNN,SUMCHN,BGAREA,LIMERR)
C
C----- USGS VERSION 3; AUG99 LAK3FM
C
C     ******************************************************************
C     ADD LAKE TERMS TO RHS AND HCOF IF SEEPAGE OCCURS IN MODEL CELLS
C     ******************************************************************
C
      COMMON /DISCOM/LBOTM(200),LAYCBD(200)
C     ------------------------------------------------------------------
C     SPECIFICATIONS:
C     ------------------------------------------------------------------
      DOUBLE PRECISION HNEW,BOTLK,BOTCL,CONDUC,H,FLOBOT,STGON,FLOBO2,
     1FLOBO1,H1,THET1
      DIMENSION CNDFCT(MXLKND),ILAKE(5,MXLKND),HNEW(NCOL,NROW,NLAY),
     1 HCOF(NCOL,NROW,NLAY),RHS(NCOL,NROW,NLAY),IBOUND(NCOL,NROW,NLAY),
     2 BOTM(NCOL,NROW,0:NBOTM),STGOLD(NLAKES)
     3 ,STRIN(NSS),STROUT(NSS),ITRB(NLAKES,NSS),IDIV(NLAKES,NSS),
     4 STGNEW(NLAKES),PRCPLK(NLAKES),EVAPLK(NLAKES),DELR(NCOL),
     5 DELC(NROW),EVAP(NLAKES),PRECIP(NLAKES),SEEP(NLAKES),SURFA(NLAKES)
     7 ,SURFIN(NLAKES),SURFOT(NLAKES),WTHDRW(NLAKES),RNF(NLAKES),
     8 DSTROT(NSS),NCNCVR(NLAKES),DSRFOT(NLAKES),SSMN(NLAKES),
     9 SUMCNN(NLAKES),SUMCHN(NLAKES),BGAREA(NLAKES),SSMX(NLAKES),
     * LIMERR(NLAKES)
C     -------------------------------------------------------------------
C1------IF LKNODE<=0 THERE ARE NO LAKE NODES. RETURN.
      IF (LKNODE.LE.0) RETURN
C
C2------PROCESS EACH CELL IN THE ILAKE LIST.
      DO 100 LK=1,NLAKES
        NCNCVR(LK) = 0
        LIMERR(LK) = 0
        SURFIN(LK)=0.0
        DSRFOT(LK)=0.0
100     SURFOT(LK)=0.0
C
C2A --- SUM UP INFLOWS FROM INFLOWING STREAM REACHES.
      DO 200 LK=1,NLAKES
         DO 200 ITRIB=1,NTRB
            INODE=ITRB(LK,ITRIB)
            IF (INODE.LE.0) GO TO 200
            SURFIN(LK)=SURFIN(LK)+STRIN(INODE)
200   CONTINUE
C2B --- SUM UP OUTFLOWS FROM OUTFLOWING STREAM REACHES.
         DO 300 LK=1,NLAKES
         DO 300 IDV=1,NDV
            INODE=IDIV(LK,IDV)
            IF (INODE.LE.0) GO TO 300
            SURFOT(LK)=SURFOT(LK)+STROUT(INODE)
            IF(ISS.NE.0) DSRFOT(LK) = DSRFOT(LK) + DSTROT(INODE)
300      CONTINUE
         ISS1 = ISS
  350    MTER = 1
         THET1 = THETA
         IF(ISS1.NE.0) MTER = NSSITR
         DO 1001 L1=1,MTER
C2C ---- INITIALIZE SUMMATION PARAMETERS.
          DO 400 LK=1,NLAKES
             SUMCNN(LK) = 0.0
             SUMCHN(LK) = 0.0
             EVAP(LK)=0.0
             PRECIP(LK)=0.0
             SEEP(LK)=0.0
400          SURFA(LK)=0.0
C
C   MASTER NODE LOOP -- COMPUTE LAKEBED SEEPAGE TERMS AND ADD TO RESIDUAL
C   AND MATRIX TERMS FOR SOLUTION.
C
      DO 900 L=1,LKNODE
         IL=ILAKE(1,L)
         IR=ILAKE(2,L)
         IC=ILAKE(3,L)
C
C4------DETERMINE LAKE AND NODAL LAYER,ROW,COLUMN NUMBER.
         LAKE=ILAKE(4,L)
         IF(ISS1.NE.0.AND.NCNCVR(LAKE).NE.0) GO TO 900
C  Convert ILAKE(5,L):  1 and 2 are type 1,  3 and 4 are type 2, 6 is type 0
         ITYPE = (ILAKE(5,L)+1)/2
         IF(ITYPE.EQ.3) ITYPE=0
         AREA = DELR(IC)*DELC(IR)
         IF(IL.GT.1) BOTLK = BOTM(IC,IR,LBOTM(IL-1))
         BOTCL = BOTM(IC,IR,LBOTM(IL))
         IF(IL.EQ.NLAY.AND.CNDFCT(L).EQ.0.0) BOTLK = BOTCL
         RAIN = PRCPLK(LAKE)
         EV = EVAPLK(LAKE)
C
C5------CONDUCTANCE FACTOR NEEDED FOR SEEPAGE CALCULATIONS.
         FLOBOT = 0.0
         CONDUC=CNDFCT(L)
         IF(CONDUC.EQ.0.0) GO TO 600
         H=HNEW(IC,IR,IL)
         STGON = (1.0-THET1)*STGOLD(LAKE) + THET1*STGNEW(LAKE)
C
C6------COMPUTE SEEPAGE INTO OR OUT OF A LAKE BED NODE, AND
C       SEEPAGE INTO OR OUT OF A LAKE WALL NODE, WHEN ITYPE=1 OR 2.
         IF (ITYPE.NE.0) GO TO 515
            IL1 = IL
            IF(IBOUND(IC,IR,IL).GT.0) GO TO 508
            DO 505 LI=IL,NLAY
            IF(IBOUND(IC,IR,LI).GT.0) GO TO 507
  505       CONTINUE
            WRITE(IOUT,506) L,IC,IR,IL
  506       FORMAT(1X,'ERROR - NO AQUIFER UNDER LAKE CELL ',4I5)
  507       IL1 = LI
            H = BOTLK
  508      IF(STGON.LE.BOTLK) GO TO 547
           IF(H.GE.BOTLK) FLOBO2=CONDUC*(STGNEW(LAKE)-H)
           IF(H.GE.BOTLK) FLOBO1=CONDUC*(STGOLD(LAKE)-H)
           IF(H.LT.BOTLK) FLOBO2=CONDUC*(STGNEW(LAKE)-BOTLK)
           IF(H.LT.BOTLK) FLOBO1=CONDUC*(STGOLD(LAKE)-BOTLK)
           FLOBOT = THET1*FLOBO2 + (1.0-THET1)*FLOBO1
           GO TO 535
  515     IF (ITYPE.NE.1.AND.ITYPE.NE.2) GO TO 550
          IF(IBOUND(IC,IR,IL).LE.0) GO TO 900
          IF(STGON.LE.BOTCL) GO TO 550
          HD = H
          IF(H.GT.BOTM(IC,IR,LBOTM(IL)-1))
     1                         HD = BOTM(IC,IR,LBOTM(IL)-1)
          THCK = HD - BOTCL
          IF(THCK.LE.0.0) THCK = 0.0
          CONDUC = CONDUC*THCK
          FLOBO2 = CONDUC*(STGNEW(LAKE)-H)
          FLOBO1 = CONDUC*(STGOLD(LAKE)-H)
          FLOBOT = THET1*FLOBO2 + (1.0-THET1)*FLOBO1
  535     SUMCNN(LAKE) = SUMCNN(LAKE) + CONDUC
          H1 = H
          IF(ITYPE.EQ.0.AND.BOTLK.GT.H) H1 = BOTLK
          SUMCHN(LAKE) = SUMCHN(LAKE) + CONDUC*H1
          GO TO 550
  547     IBOUND(IC,IR,IL-1) = 0
  550     CONTINUE
C
C  ADD LAKE SEEPAGE RATES TO MATRIX AND RESIDUAL TERMS BEFORE SOLVING
C
       IF(ISS1.NE.0) GO TO 600
         IF (ITYPE.NE.0) GO TO 825
           IF (STGON.LE.BOTLK) GO TO 600
           IF (H.LE.BOTLK) GO TO 870
           IF (H.GT.BOTLK) GO TO 850
825      IF (ITYPE.NE.1.AND.ITYPE.NE.2) GO TO 600
            IF (STGON.GT.BOTCL.AND.THCK.GT.0.0) GO TO 850
            GO TO 600
850      RHS(IC,IR,IL)=RHS(IC,IR,IL) - STGON*CONDUC
         HCOF(IC,IR,IL)=HCOF(IC,IR,IL) - CONDUC
         GO TO 600
870      RHS(IC,IR,IL1)=RHS(IC,IR,IL1)-FLOBOT
600      CONTINUE
C
C7------COMPUTE EVAPORATION FOR EACH LAKE NODE AND ADD IT TO THE
C          ACCUMULATIVE EVAPORATION DISCHARGE.
          IF (ITYPE.EQ.0) THEN
             EP=EV*AREA
             IF (STGON.LE.BOTLK) EP=0.0
             EVAP(LAKE)=EVAP(LAKE)+EP
          ENDIF
C8-------COMPUTE RAINFALL RECHARGE FOR EACH LAKE NODE AND ADD IT TO
C         THE ACCUMULATIVE PRECIPITATION RECHARGE.
          IF (ITYPE.EQ.0) THEN
             R2 = RAIN*AREA
             IF(STGON.LE.BOTLK) R2 = 0.0
             PRECIP(LAKE)=PRECIP(LAKE)+R2
          ENDIF
C9-------COMPUTE ACCUMULATIVE SEEPAGE THROUGH LAKEBED SEDIMENTS.
          IF(ISS1.NE.0) SEEP(LAKE)=SEEP(LAKE)-FLOBOT
C10------COMPUTE ACCUMULATIVE LAKE NODAL AREA, EXCLUDE WALL NODES.
          IF (ITYPE.EQ.0.AND.STGON.GT.BOTLK) THEN
              SURFA(LAKE)=SURFA(LAKE)+AREA
          ENDIF
900      CONTINUE
C
C11------ONLY COMPUTE LAKE LEVEL AFTER SCANNING THRU ALL NODES OF A LAKE.
C
       DO 1000 LAKE=1,NLAKES
       IF(SURFA(LAKE).LE.0.0) GO TO 1000
       IF(RNF(LAKE).GE.0.0) RUNF = RNF(LAKE)
       IF(RNF(LAKE).LT.0.0) RUNF =-RNF(LAKE)*PRCPLK(LAKE)*BGAREA(LAKE)
       IF(ISS1.NE.0) GO TO 970
       DDSF = DELT/SURFA(LAKE)
       STGNEW(LAKE)=(STGOLD(LAKE)+DDSF*(PRECIP(LAKE)
     1 -EVAP(LAKE)+RUNF-WTHDRW(LAKE)+SURFIN(LAKE)-SURFOT(LAKE)-
     2 (1.0-THETA)*STGOLD(LAKE)*SUMCNN(LAKE)+
     3 SUMCHN(LAKE)))/(1.0+THETA*DDSF*SUMCNN(LAKE))
      GO TO 1000
  970 IF(NCNCVR(LAKE).NE.0) GO TO 1000
      SSMN1 = SSMN(LAKE)
      SSMX1 = SSMX(LAKE)
      RESID = PRECIP(LAKE)-EVAP(LAKE)+RUNF-WTHDRW(LAKE)+
     1  SURFIN(LAKE)-SURFOT(LAKE)+SEEP(LAKE)
      STGO = STGOLD(LAKE)
      STGOLD(LAKE) = STGOLD(LAKE) + RESID/(DSRFOT(LAKE)+SUMCNN(LAKE))
      IF(STGOLD(LAKE).LE.SSMX1.AND.STGOLD(LAKE).GE.SSMN1) GO TO 975
      WRITE(IOUT,976) STGOLD(LAKE)
  976 FORMAT(/1X,'ERROR -- COMPUTED STEADY-STATE STAGE EXCEEDS SPECIFIED
     1 MAXIMUM OR MINIMUM ',F10.2)
      LIMERR(LAKE) = LIMERR(LAKE) + 1
      IF(LIMERR(LAKE).LT.10) GO TO 977
      WRITE(IOUT,979) LAKE
  979 FORMAT(1X,'LAKE',I3,3X,'STAGE REPEATEDLY EXCEEDS ITS BOUNDS -- STO
     1PPING EXECUTION')
        call stopfile ! emrl
      STOP
  977 STGOLD(LAKE) = (SSMN1+SSMX1)/2.0
  975 STGNEW(LAKE) = STGOLD(LAKE)
      DSTG = ABS(STGOLD(LAKE)-STGO)
      IF(DSTG.LE.SSCNCR) NCNCVR(LAKE) = 1
 1000  CONTINUE
       IF(ISS1.EQ.0) RETURN
       NCNV = 0
       DO 1002 LAKE=1,NLAKES
       IF(NCNCVR(LAKE).EQ.0.AND.SURFA(LAKE).GT.0.0) NCNV = 1
 1002  CONTINUE
       IF(NCNV.EQ.0) GO TO 1005
 1001  CONTINUE
       DO 1003 LAKE=1,NLAKES
       IF(NCNCVR(LAKE).EQ.0) WRITE(IOUT,1004) KITER, LAKE, STGOLD(LAKE)
 1004  FORMAT(1X,'ITERATION ',I3,2X,'LAKE ',I3,2X,'STAGE ',F10.4,2X,
     1  'DID NOT CONVERGE TO STEADY-STATE')
 1003  CONTINUE
       GO TO 350
 1005  WRITE(IOUT,1006) L1
 1006  FORMAT(1X,'ALL LAKES CONVERGED TO STEADY-STATE AFTER ',I3,2X,
     1  'ITERATIONS')
       WRITE(IOUT,1007) (STGNEW(LAKE),LAKE=1,NLAKES)
 1007  FORMAT(1X,'LAKE STAGES = ',10F11.3)
       ISS1 = 0
       GO TO 350
C
      END
C
      SUBROUTINE LAK3BD(LKNODE,MXLKND,NODES,ILAKE,HNEW,
     1      IBOUND,NCOL,NROW,NLAY,NLAKES,DELT,NSS,NTRB,NDV,
     2      ITRB,IDIV,STRIN,STROUT,STGOLD,STGNEW,CNDFCT,
     3      PRCPLK,EVAPLK,DELR,DELC,BOTM,NBOTM,VBVL,VBNM,MSUM,
     4      KSTP,KPER,ILKCB,ICBCFL,BUFF,IOUT,STAGES,PERTIM,TOTIM,
     5      ICS,ISUB,SILLVT,ICMX,NCLS,WTHDRW,LWRT,
     6      EVAP,PRECIP,SEEP,SURFA,SURFIN,SURFOT,VOL,GWIN,GWOUT,DELH,
     7      TDELH,LDRY,BOTTMS,NCNT,KSUB,STGADJ,FLXINL,NCNST,SVT,JCLS,
     8      RNF,THETA,SUMCNN,SUMCHN,BGAREA,IUNITGWT,KCNT,MSUB,MSUB1,
     9      IUNITGAGE,NUMGAGE,IGGLST,VOLOLD,FLOB,ISS,LAYHDT,
     9      IAUXSV,VOLINIT,STGOLD2,CUMPPT,CUMEVP,CUMRNF,
     1      CUMGWI,CUMGWO,CUMSWI,CUMSWO,CUMWDR,CUMFLX)
C
C----- USGS VERSION 3; AUG 1999 LAK3BD
C
C     ******************************************************************
C     CALCULATE VOLUMETRIC BUDGET FOR LAKES
C     ******************************************************************
C
C     ------------------------------------------------------------------
C     SPECIFICATIONS:
C     ------------------------------------------------------------------
      COMMON /DISCOM/LBOTM(200),LAYCBD(200)
      CHARACTER*16 VBNM,TEXT
      DOUBLE PRECISION HNEW,BOTLK,BOTCL,CONDUC,H,FLOBOT,STGON,FLOBO2,
     1FLOBO1,H1,RATE,RATIN,RATOUT
      DIMENSION CNDFCT(MXLKND),ILAKE(5,MXLKND),HNEW(NCOL,NROW,NLAY),
     1 IBOUND(NCOL,NROW,NLAY),RNF(NLAKES),BGAREA(NLAKES),
     2 STRIN(NSS),STROUT(NSS),ITRB(NLAKES,NSS),IDIV(NLAKES,NSS),
     3 STGOLD(NLAKES),STGNEW(NLAKES),PRCPLK(NLAKES),EVAPLK(NLAKES),
     4 DELR(NCOL),DELC(NROW),BOTM(NCOL,NROW,0:NBOTM),LAYHDT(NLAY),
     5 VBVL(4,MSUM),VBNM(MSUM),BUFF(NCOL,NROW,NLAY),STAGES(NLAKES),
     6 ICS(NLAKES),ISUB(NLAKES,NLAKES),SILLVT(NLAKES,NLAKES)
       DIMENSION JCLS(NCLS,ICMX),NCNT(NLAKES),KSUB(NLAKES),
     7 STGADJ(NLAKES),MSUB(NLAKES,NLAKES),MSUB1(NLAKES),
     8 FLXINL(NLAKES),NCNST(NLAKES),SVT(NLAKES),WTHDRW(NLAKES),
     9 IGGLST(3,NUMGAGE),VOLOLD(NLAKES)
      DIMENSION FLOB(LKNODE)
      DIMENSION EVAP(NLAKES),PRECIP(NLAKES),SEEP(NLAKES),
     1 SURFA(NLAKES),SURFIN(NLAKES),SURFOT(NLAKES),VOL(NLAKES),
     2 GWIN(NLAKES),GWOUT(NLAKES),DELH(NLAKES),TDELH(NLAKES),
     3 LDRY(NODES),BOTTMS(NLAKES),SUMCNN(NLAKES),SUMCHN(NLAKES)
      DIMENSION CUMPPT(NLAKES),CUMEVP(NLAKES),CUMRNF(NLAKES),
     1 CUMGWI(NLAKES),CUMGWO(NLAKES),CUMSWI(NLAKES),CUMSWO(NLAKES),
     2 CUMWDR(NLAKES),CUMFLX(NLAKES)
cgage
      DIMENSION VOLINIT(NLAKES),STGOLD2(NLAKES)
      DIMENSION ILB(5),IRB(5),ICB(5)
      CHARACTER*16 LAKAUX(5)
      DIMENSION FACE(1)
      DATA TEXT /'   LAKE  SEEPAGE'/
      DATA LAKAUX(1)/'IFACE'/
C     ------------------------------------------------------------------
C
C1------SET IBD=1 IF BUDGET TERMS SHOULD BE SAVED ON DISK.
      ZERO = 0.0
      IBD=0
      KCNT = 0
      RATIN = 0.
      RATOUT =0.
      DO 104 LDR = 1,NODES
  104 LDRY(LDR) = 0
      LDR = 0
C
C1A-----TEST TO SEE IF CELL-BY-CELL TERMS ARE NEEDED.
      IF(ILKCB.GT.0) IBD=ICBCFL
C  -----IF COMPACT BUDGET, WRITE LIST HEADER
      IF(IBD.EQ.2) THEN
         NAUX=0
         IF(IAUXSV.NE.0) NAUX=1
         CALL UBDSV4(KSTP,KPER,TEXT,NAUX,LAKAUX,ILKCB,NCOL,NROW,NLAY,
     1               LKNODE,IOUT,DELT,PERTIM,TOTIM,IBOUND)
      END IF
C
C1B------IF NO LAKE NODES, KEEP ZERO IN ACCUMULATORS.
      IF (LKNODE.EQ.0) GO TO 1200
C
C1C-----CLEAR CELL-BY-CELL BUFFER.
      DO 5 IL=1,NLAY
      DO 5 IR=1,NROW
      DO 5 IC=1,NCOL
5        BUFF(IC,IR,IL)=ZERO
C------
C2------PROCESS EACH CELL IN THE ILAKE LIST.
      DO 100 LK=1,NLAKES
        FLXINL(LK)=ZERO
        SURFIN(LK)=ZERO
100     SURFOT(LK)=ZERO
C
CMOC---If transport active, save previous value of lake volume.
      IF (IUNITGWT.GT.0.OR.NUMGAGE.GT.0) THEN
        DO 150 LK=1,NLAKES
          VOLOLD(LK)=VOL(LK)
150     CONTINUE
      END IF
C
C2A --- SUM UP INFLOWS FROM INFLOWING STREAM REACHES.
      DO 200 LK=1,NLAKES
         DO 200 ITRIB=1,NTRB
            INODE=ITRB(LK,ITRIB)
            IF (INODE.LE.0) GO TO 200
            SURFIN(LK)=SURFIN(LK)+STRIN(INODE)*DELT
200   CONTINUE
C2B --- SUM UP OUTFLOWS FROM OUTFLOWING STREAM REACHES.
         DO 300 LK=1,NLAKES
         DO 300 IDV=1,NDV
            INODE=IDIV(LK,IDV)
            IF (INODE.LE.0) GO TO 300
            SURFOT(LK)=SURFOT(LK)+STROUT(INODE)*DELT
300      CONTINUE
C2C ---- INITIALIZE SUMMATION PARAMETERS.
          DO 400 LK=1,NLAKES
             SUMCNN(LK) = ZERO
             SUMCHN(LK) = ZERO
             EVAP(LK)=ZERO
             PRECIP(LK)=ZERO
             SEEP(LK)=ZERO
             VOL(LK)=ZERO
400          SURFA(LK)=ZERO
        DO 425 NM=1,NLAKES
           GWIN(NM)=ZERO
425       GWOUT(NM)=ZERO
C
C   MASTER NODE LOOP -- COMPUTE LAKEBED SEEPAGE TERMS AND BUDGET TERMS
C
         IF (ILKCB.LT.0.AND.ICBCFL.NE.0) WRITE (IOUT,'(//)')
      DO 900 L=1,LKNODE
         IL=ILAKE(1,L)
         IR=ILAKE(2,L)
         IC=ILAKE(3,L)
         RATE=ZERO
C
C4------DETERMINE LAKE AND NODAL LAYER,ROW,COLUMN NUMBER.
         LAKE=ILAKE(4,L)
C  Convert ILAKE(5,L):  1 and 2 are type 1,  3 and 4 are type 2, 6 is type 0
         ITYPE = (ILAKE(5,L)+1)/2
         IF(ITYPE.EQ.3) ITYPE=0
         AREA = DELR(IC)*DELC(IR)
         IF(IL.GT.1) BOTLK = BOTM(IC,IR,LBOTM(IL-1))
         BOTCL = BOTM(IC,IR,LBOTM(IL))
         IF(IL.EQ.NLAY.AND.CNDFCT(L).EQ.ZERO) BOTLK = BOTCL
         RAIN = PRCPLK(LAKE)
         EV = EVAPLK(LAKE)
C5------CONDUCTANCE FACTOR NEEDED FOR SEEPAGE CALCULATIONS.
         FLOBOT = ZERO
         CONDUC=CNDFCT(L)
         IF(CONDUC.EQ.ZERO) GO TO 550
         H=HNEW(IC,IR,IL)
         STGON = (1.0-THETA)*STGOLD(LAKE) + THETA*STGNEW(LAKE)
C
C6------COMPUTE SEEPAGE INTO OR OUT OF A LAKE BED NODE, AND
C       SEEPAGE INTO OR OUT OF A LAKE WALL NODE, WHEN ITYPE=1 OR 2.
         IF (ITYPE.NE.0) GO TO 515
            IF(IBOUND(IC,IR,IL).GT.0) GO TO 508
            DO 505 LI=IL,NLAY
            IF(IBOUND(IC,IR,LI).GT.0) GO TO 507
  505       CONTINUE
            WRITE(IOUT,506) L,IC,IR,IL
  506       FORMAT(1X,'ERROR - NO AQUIFER UNDER LAKE CELL ',4I5)
        call stopfile ! emrl
            STOP
  507       H = BOTLK
  508      IF(STGON.LE.BOTLK) GO TO 550
           IF(H.GE.BOTLK) FLOBO2=CONDUC*(STGNEW(LAKE)-H)
           IF(H.GE.BOTLK) FLOBO1=CONDUC*(STGOLD(LAKE)-H)
           IF(H.LT.BOTLK) FLOBO2=CONDUC*(STGNEW(LAKE)-BOTLK)
           IF(H.LT.BOTLK) FLOBO1=CONDUC*(STGOLD(LAKE)-BOTLK)
           FLOBOT = THETA*FLOBO2 + (1.0-THETA)*FLOBO1
           GO TO 535
  515     IF (ITYPE.NE.1.AND.ITYPE.NE.2) GO TO 550
          IF(IBOUND(IC,IR,IL).LE.0) GO TO 899
          IF(STGON.LE.BOTCL) GO TO 550
          HD = H
          IF(H.GT.BOTM(IC,IR,LBOTM(IL)-1))
     1                      HD = BOTM(IC,IR,LBOTM(IL)-1)
          THCK = HD - BOTCL
          IF(THCK.LE.ZERO) THCK = ZERO
          CONDUC = CONDUC*THCK
          FLOBO2 = CONDUC*(STGNEW(LAKE)-H)
          FLOBO1 = CONDUC*(STGOLD(LAKE)-H)
          FLOBOT = THETA*FLOBO2 + (1.0-THETA)*FLOBO1
  535     SUMCNN(LAKE) = SUMCNN(LAKE) + CONDUC
          H1 = H
          IF(ITYPE.EQ.0.AND.BOTLK.GT.H) H1 = BOTLK
          SUMCHN(LAKE) = SUMCHN(LAKE) + CONDUC*H1
  550     CONTINUE
C
C7------COMPUTE EVAPORATION FOR EACH LAKE NODE AND ADD IT TO THE
C          ACCUMULATIVE EVAPORATION DISCHARGE.
          IF (ITYPE.EQ.0) THEN
             EP=EV*AREA
             IF (STGON.LE.BOTLK) EP=ZERO
             EVAP(LAKE)=EVAP(LAKE)+EP*DELT
          ENDIF
C
C8-------COMPUTE RAINFALL RECHARGE FOR EACH LAKE NODE AND ADD IT TO
C         THE ACCUMULATIVE PRECIPITATION RECHARGE.
          IF (ITYPE.EQ.0) THEN
             R2 = RAIN*AREA
             IF(STGON.LE.BOTLK) R2 = ZERO
             PRECIP(LAKE)=PRECIP(LAKE)+R2*DELT
          ENDIF
C
C9-------COMPUTE ACCUMULATIVE SEEPAGE THROUGH LAKEBED SEDIMENTS.
C         SEEP(LAKE)=SEEP(LAKE)-FLOBOT*DELT
C
C10------COMPUTE ACCUMULATIVE LAKE NODAL AREA, EXCLUDE WALL NODES.
          IF (ITYPE.EQ.0.AND.STGON.GT.BOTLK) THEN
             SURFA(LAKE)=SURFA(LAKE)+AREA
             VOL(LAKE)=VOL(LAKE)+(STGON-BOTLK)*AREA
          ENDIF
          RATE=FLOBOT
          IF (IUNITGWT.GT.0) FLOB(L)=FLOBOT
          IF (ILKCB.LT.0.AND.ICBCFL.NE.0) WRITE(IOUT,880)
     1          TEXT,KPER,KSTP,L,IL,IR,IC,RATE
880       FORMAT(1X,A,'   PERIOD',I3,'   STEP',I3,'   NODE',I4,
     1            '   LAYER',I3,'   ROW',I4,'   COL',I4,'   RATE',
     2            G15.7)
C
C------ ADD RATE TO BUFFER.
          BUFF(IC,IR,IL)=BUFF(IC,IR,IL)+RATE
C
C------ SEE IF FLOW IS INTO AQUIFER OR INTO LAKE.
            IF (RATE) 885,899,890
C
C------ AQUIFER IS DISCHARGING TO LAKE SUBTRACT RATE FROM RATOUT.
885         RATOUT=RATOUT-RATE
            GWIN(LAKE)=GWIN(LAKE)-RATE
            GO TO 899
C
C------ AQUIFER IS RECHARGED FROM LAKE ADD RATE TO RATIN.
890         RATIN=RATIN+RATE
            GWOUT(LAKE)=GWOUT(LAKE)+RATE

C------ IF SAVING COMPACT BUDGET, WRITE FLOW FOR ONE LAKE FACE
899      IF(IBD.EQ.2) THEN
            FACE(1)=ILAKE(5,L)
            R=RATE
            CALL UBDSVB(ILKCB,NCOL,NROW,IC,IR,IL,R,FACE(1),1,NAUX,
     1                 1,IBOUND,NLAY)
         END IF
900      CONTINUE
C
C11------ONLY COMPUTE LAKE LEVEL AFTER SCANNING THRU ALL NODES OF A LAKE.
C
      IF(ISS.LE.0) GO TO 905
      DO 903 LAKE=1,NLAKES
            DELH(LAKE)=STGNEW(LAKE)-STGOLD(LAKE)
            TDELH(LAKE)=STGNEW(LAKE)-STAGES(LAKE)
  903 CONTINUE
      GO TO 1350
  905     DO 1000 LAKE=1,NLAKES
         STGON = (1.0-THETA)*STGOLD(LAKE) + THETA*STGNEW(LAKE)
          WDRAW=WTHDRW(LAKE)*DELT
       IF(RNF(LAKE).GE.ZERO) RUNF = RNF(LAKE)
       IF(RNF(LAKE).LT.ZERO) RUNF =-RNF(LAKE)*PRCPLK(LAKE)*BGAREA(LAKE)
         RUNFD = RUNF*DELT
          IF(VOL(LAKE).LE.0) GO TO 1110
       STGNEW(LAKE)=(STGOLD(LAKE)*SURFA(LAKE)+(PRECIP(LAKE)-EVAP(LAKE)
     1  -WDRAW+RUNFD+SURFIN(LAKE)-SURFOT(LAKE)-
     2 (1.0-THETA)*DELT*STGOLD(LAKE)*SUMCNN(LAKE)+
     3 DELT*SUMCHN(LAKE)))/(SURFA(LAKE)+THETA*DELT*SUMCNN(LAKE))
         VOL(LAKE)=VOL(LAKE) + (STGNEW(LAKE)-STGON)*SURFA(LAKE)
      IF(VOL(LAKE).LE.ZERO) WRITE(IOUT,1114) LAKE
 1114 FORMAT(1X,'..........LAKE',I3,' HAS JUST GONE DRY..........')
            DELH(LAKE)=STGNEW(LAKE)-STGOLD(LAKE)
            TDELH(LAKE)=STGNEW(LAKE)-STAGES(LAKE)
      GO TO 1000
C
C------CHECK FOR CONDITION (AVERAGE HEAD IN UNDERLYING AQUIFER
C      GREATER THAN BOTTOM OF LAKE) FOR REWETTING A DRY LAKE
C
 1110 AVHD = ZERO
      BOTARE = ZERO
      WRITE(IOUT,1112) LAKE
 1112 FORMAT(1X,'..........LAKE',I3,' IS DRY..........')
      IF(NLAY.EQ.1) GO TO 1000
      DO 1115 L=1,LKNODE
      L1 = ILAKE(4,L)
C  Convert ILAKE(5,L):  1 and 2 are type 1,  3 and 4 are type 2, 6 is type 0
      ITYPE = (ILAKE(5,L)+1)/2
      IF(ITYPE.EQ.3) ITYPE=0
      IF(L1.NE.LAKE.OR.ITYPE.NE.0) GO TO 1115
      K = ILAKE(1,L)
      J = ILAKE(2,L)
      I = ILAKE(3,L)
      IF(K.EQ.NLAY.AND.IBOUND(I,J,K).EQ.0) GO TO 1000
      IF(K.EQ.1) GO TO 1115
      IF(BOTM(I,J,LBOTM(K-1)).GT.BOTTMS(LAKE)) GO TO 1115
      DO 1117 K1=K,NLAY
      IF(LAYHDT(K1).EQ.0) GO TO 1117
      K2 = K1
      IF(IBOUND(I,J,K1).LE.0) GO TO 1117
      GO TO 1119
 1117 CONTINUE
      AVHD = AVHD + BOTM(I,J,LBOTM(K2))*DELR(I)*DELC(J)
      GO TO 1121
 1119 AVHD = AVHD + HNEW(I,J,K1)*DELR(I)*DELC(J)
 1121 BOTARE = BOTARE + DELR(I)*DELC(J)
 1115 CONTINUE
      IF(BOTARE.LE.ZERO) GO TO 1000
      AVHD = AVHD/BOTARE
      IF(AVHD.LT.BOTTMS(LAKE)) GO TO 1125
      WRITE(IOUT,1122)
 1122 FORMAT(/1X,'AQUIFER HEAD UNDERNEATH BOTTOM OF LAKE IS HIGHER THAN
     1LAKE BOTTOM ELEVATION')
 1120 STGNEW(LAKE) = AVHD
      SURFA(LAKE) = BOTARE
      VOL(LAKE) = (STGNEW(LAKE)-BOTTMS(LAKE))*SURFA(LAKE)
      WRITE(IOUT,1184) LAKE, STGNEW(LAKE)
 1184 FORMAT(/1X,'LAKE',I3,' HAS REWET.  SET STAGE OF LAKE TO',F10.2,
     1  2X,'FT'/)
      GO TO 1000
C
C-----CHECK FOR STREAM OR AUGMENTATION INFLOWS
C
 1125 TSTVOL = SURFIN(LAKE)
      IF(WDRAW.LT.ZERO) TSTVOL = TSTVOL - WDRAW
      IF(TSTVOL.LE.ZERO) GO TO 1000
      STGNEW(LAKE) = BOTTMS(LAKE)
      SURFA(LAKE) = BOTARE
      AVHD = BOTTMS(LAKE) + TSTVOL/BOTARE
      WRITE(IOUT,1123)
 1123 FORMAT(/1X,'THERE ARE NON-ZERO SURFACE-WATER INFLUXES INTO THE DRY
     1 LAKE')
      GO TO 1120
C
 1000     CONTINUE
C
C-----ADJUST STAGES OF COALESCENT MULTIPLE-LAKE SYSTEMS
C
      KCNT = 0
      IF(NCLS.LE.0) GO TO 1350
      DO 1205 I=1,NCLS
      DO 1205 J=1,ICMX
 1205 JCLS(I,J) = 0
C
C   CHECK EACH LAKE SYSTEM (ICL) FOR CURRENT CONNECTIONS TO SUBLAKES
C
      DO 1300 ICL=1,NCLS
      DO 1206 K=1,NLAKES
      SVT(K) = ZERO
      NCNST(K) = 0
 1206 NCNT(K) = 0
C
C   ELIMINATE (JCLS=2) ALL LAKES THAT HAVE ALREADY HAD THEIR STAGES
C   ADJUSTED AS PART OF A CONNECTED SYSTEM OF LAKES AND SUBLAKES
C
      DO 1210 IC4=ICL,NCLS
      ICM4 = ICS(IC4)
      DO 1210 IC5=1,ICM4
      IF(JCLS(IC4,IC5).EQ.1) JCLS(IC4,IC5) = 2
 1210 CONTINUE
 1215 IF(JCLS(ICL,1).GE.2) GO TO 1300
C
C   TAG CENTER LAKE BY SETTING JCLS=1 AND THEN CHECK SUBLAKES FOR
C   CONNECTIONS.  IF CONNECTED, SET JCLS=1 AND NCNT=1.
C
      ICM = ICS(ICL)
      IS1 = ISUB(ICL,1)
      JCLS(ICL,1) = 1
      NCNT(IS1) = 1
      SVT(IS1) = -99999.
      DO 1220 J=2,ICM
      IS2 = ISUB(ICL,J)
      IF(IS2.LE.0) GO TO 1225
      IF(STGNEW(IS1).LE.SILLVT(ICL,J-1).AND.STGNEW(IS2).LE.
     1  SILLVT(ICL,J-1)) GO TO 1220
      JCLS(ICL,J) = 1
      NCNT(IS2) = 1
      SVT(IS2) = SILLVT(ICL,J-1)
 1220 CONTINUE
 1225 IF(ICL.EQ.NCLS) GO TO 1240
C
C   CHECK TO SEE IF CENTER LAKES OF REMAINING LAKE SYSTEMS ARE THE SAME
C   AS CONNECTED SUBLAKES OF THE PRESENT LAKE SYSTEM.  IF SO, CHECK
C   THEIR SUBLAKES FOR CONNECTIONS TO THE CENTER LAKES.  IF SUBLAKES
C   ARE ADDED TO THE SYSTEM, THEN CHECK REMAINING CENTER LAKES FOR AN
C   IDENTITY WITH THE NEWLY ADDED SUBLAKES.  ALL CONNECTED LAKES ARE
C   APPROPRIATELY TAGGED (JCLS=1 AND NCNT=1).
C
      ICL1 = ICL + 1
      DO 1230 LK=ICL1,NCLS
      IF(JCLS(LK,1).EQ.2) GO TO 1230
      LK3 = LK - 1
      DO 1227 LK1=1,LK3
      ICM2 = ICS(LK1)
      DO 1227 IC2=2,ICM2
      IF(ISUB(LK1,IC2).NE.ISUB(LK,1).OR.JCLS(LK1,IC2).NE.1) GO TO 1227
      JCLS(LK,1) = 1
      IS1 = ISUB(LK,1)
      NCNT(IS1) = 1
      ICM3 = ICS(LK)
      DO 1226 IC3=2,ICM3
      IS2 = ISUB(LK,IC3)
      IF(IS2.LE.0) GO TO 1227
      IF(STGNEW(IS1).LE.SILLVT(LK,IC3-1).AND.STGNEW(IS2).LE.
     1  SILLVT(LK,IC3-1)) GO TO 1226
      JCLS(LK,IC3) = 1
      NCNT(IS2) = 1
      SVT(IS2) = SILLVT(LK,IC3-1)
 1226 CONTINUE
 1227 CONTINUE
 1230 CONTINUE
C
C   COUNT NUMBER OF LAKES IDENTIFIED AS A CONNECTED PART OF THE PRESENT
C   LAKE SYSTEM, STORE LAKE NUMBERS IN KSUB, AND CUMULATE TOTAL SURFACE
C   AREA
C
 1240 ICNT = 0
      KCNT = KCNT + 1
      TOTARE = ZERO
      DO 1245 L=1,NLAKES
      IF(NCNT(L).NE.1) GO TO 1245
      ICNT = ICNT + 1
      TOTARE = TOTARE + SURFA(L)
      KSUB(ICNT) = L
      MSUB(ICNT,KCNT) = L
 1245 CONTINUE
      MSUB1(KCNT) = ICNT
      IF(ICNT.LE.1) KCNT = KCNT - 1
      IF(ICNT.LE.1) GO TO 1300
      IF(LWRT.GT.0.OR.ICBCFL.LE.0) GO TO 1251
      WRITE(IOUT,1250) KSTP, ICNT, TOTARE
 1250 FORMAT(/1X,80('-')/1X,'TIME STEP ',I3,5X,'NUMBER OF CONNECTED LAKE
     1S IS',I3,5X,'TOTAL AREA = ',D16.9/)
 1251 CONTINUE
C
C   COMPUTE STAGE ADJUSTMENTS (STGADJ) REQUIRED FOR CONNECTED LAKES TO
C   HAVE THE SAME STAGE.
C
      DO 1270 I=1,ICNT
      L1 = KSUB(I)
      SUM = ZERO
      DO 1265 J=1,ICNT
      IF(J.EQ.I) GO TO 1265
      L = KSUB(J)
      SUM = SUM + SURFA(L)*(STGNEW(L)-STGNEW(L1))
 1265 CONTINUE
      STGADJ(L1) = SUM/TOTARE
 1270 CONTINUE
C
C   CHECK FOR NEWLY COMPUTED LAKE STAGES (STGNEW) LESS THAN SILL
C   ELEVATIONS OF THE LAKES (SVT)
C
 1272 ICNR = 0
      ICM3 = ICS(ICL)
      DO 1275 IC3=2,ICM3
      L = ISUB(ICL,IC3)
      IF(L.LE.0) GO TO 1275
      IF(NCNST(L).EQ.1) GO TO 1275
      IF(NCNT(L).EQ.0) GO TO 1275
      IF(STGNEW(L).LT.SVT(L)) GO TO 1275
      STGTST = STGNEW(L) + STGADJ(L)
      IF(STGTST.GE.SVT(L)) GO TO 1275
C
C   ADJUST STAGE TO SILL ELEVATION
C
      NCNST(L) = 1
      STGADJ(L) = SVT(L) - STGNEW(L)
      STGNEW(L) = SVT(L)
      FLXINL(L) = STGADJ(L)*SURFA(L)
      VOL(L) = VOL(L) + FLXINL(L)
      TOTARE = TOTARE - SURFA(L)
      NCNT(L) = 2
      ICNR = 1
      WRITE(IOUT,2238) L,L,STGADJ(L),STGNEW(L)
 1275 CONTINUE
      IF(ICL.EQ.NCLS) GO TO 1277
C
C   IF A LAKE STAGE IS ADJUSTED TO THE SILL ELEVATION, CHECK TO SEE
C   WHETHER THERE ARE SUBLAKES OF THIS LAKE AND ADJUST THEM TO THE SILL
C   ELEVATION UNLESS THE ORIGINAL STAGE IS ALREADY LOWER, IN WHICH CASE
C   THEY ARE NO LONGER CONNECTED.
C
      ICL1 = ICL + 1
      DO 2230 LK=ICL1,NCLS
      IS1 = ISUB(LK,1)
      IF(NCNT(IS1).EQ.0) GO TO 2230
      ICM1 = ICS(LK)
      DO 2225 IC2=2,ICM1
      IS2 = ISUB(LK,IC2)
      IF(NCNST(IS2).EQ.1) GO TO 2225
      IF(NCNT(IS2).EQ.0) GO TO 2225
      IF(STGNEW(IS2).LT.SVT(IS2)) GO TO 2225
      SVT1 = SVT(IS2)
      IF(SVT(IS1).GT.SVT1.AND.NCNT(IS1).EQ.2) SVT1 = SVT(IS1)
      STGTST = STGNEW(IS2) + STGADJ(IS2)
      IF(STGTST.GE.SVT1) GO TO 2225
      ICNR = 1
      NCNST(IS2) = 1
      NCNT(IS2) = 2
      STGADJ(IS2) = SVT1 - STGNEW(IS2)
      STGNEW(IS2) = SVT1
      FLXINL(IS2) = STGADJ(IS2)*SURFA(IS2)
      VOL(IS2) = VOL(IS2) + FLXINL(IS2)
      TOTARE = TOTARE - SURFA(IS2)
      IF(LWRT.GT.0.OR.ICBCFL.LE.0) GO TO 2225
      L11 = L
      IF(SVT(IS2).GT.SVT(L)) L11 = IS2
      WRITE(IOUT,2238) IS2,L11,STGADJ(IS2),STGNEW(IS2)
 2238 FORMAT(1X,'READJUST STAGE OF LAKE ',I3,' TO LAKE ',I3,
     1  ' SILL ELEVATION BY ',F5.2,' TO ',F7.2)
 2225 CONTINUE
 2230 CONTINUE
 1277 IF(ICNR.LE.0) GO TO 1280
C
C   RECOMPUTE STAGE ADJUSTMENTS CONSTRAINED NOT TO LOWER LAKES BELOW
C   SILL ELEVATIONS.
C
      ICNR1 = 0
      DO 1370 I=1,ICNT
      L1 = KSUB(I)
      IF(NCNT(L1).EQ.0) GO TO 1370
      IF(NCNST(L1).EQ.1) GO TO 1370
      SUM = ZERO
      DO 1365 J=1,ICNT
      IF(J.EQ.I) GO TO 1365
      L = KSUB(J)
      IF(NCNT(L).EQ.0) GO TO 1365
      IF(NCNST(L).EQ.0) SUM = SUM + SURFA(L)*(STGNEW(L)-STGNEW(L1))
      IF(NCNST(L).EQ.1) SUM = SUM - SURFA(L)*STGADJ(L)
 1365 CONTINUE
      STGADJ(L1) = SUM/TOTARE
      STGTST = STGNEW(L1) + STGADJ(L1)
      IF(STGTST.GE.SVT(L1)) GO TO 1370
      ICNR1 = 1
 1370 CONTINUE
      IF(ICNR1.NE.0) GO TO 1272
 1280 IF(LWRT.GT.0.OR.ICBCFL.LE.0) GO TO 1281
      WRITE(IOUT,1286)
 1286 FORMAT(//11X,'SURFACE',7X,'SILL',3X,'WATER BUDGET',2X,'STAGE',2X,
     1 'CORRECTED',3X,'LAKE VOLUME'/2X,'LAKE',7X,'AREA',6X,'ELEVATION',
     2  3X,'STAGE',3X,'CORRECTION',2X,'STAGE',5X,'CORRECTION')
 1281 CONTINUE
      TVOLM = ZERO
      DO 1290 I=1,ICNT
      L = KSUB(I)
      STO = STGNEW(L)
      IF(NCNST(L).EQ.1) STO = STGNEW(L) - STGADJ(L)
      IF(NCNST(L).EQ.1) GO TO 1285
      STGNEW(L) = STGNEW(L) + STGADJ(L)
      FLXINL(L) = STGADJ(L)*SURFA(L)
      VOL(L) = VOL(L) + FLXINL(L)
 1327 FORMAT(/10X,'WARNING -- SUM OF INTERLAKE FLUXES ',F10.0,' EXCEEDS
     110**6 OF THE TOTAL VOLUME'/)
      WRITE(IOUT,1301)
 1301 FORMAT(1X,80('-')/)
 1285 TVOLM = TVOLM + VOL(L)
      IF(LWRT.GT.0.OR.ICBCFL.LE.0) GO TO 1290
      WRITE(IOUT,1269) L,SURFA(L),SVT(L),STO,STGADJ(L),STGNEW(L),
     1  FLXINL(L)
 1269 FORMAT(1X,I3,1X,G15.5,4F10.2,G15.5)
 1290 CONTINUE
C
C   RECOMPUTE TIME STEP AND CUMULATIVE STAGE CHANGES FOR CONNECTED LAKES
C
      DO 1295 I=1,ICNT
      L = KSUB(I)
      DELH(L) = STGNEW(L) - STGOLD(L)
      TDELH(L) = STGNEW(L) - STAGES(L)
 1295 CONTINUE
 1300 CONTINUE
C
C   CHECK ON SUM OF CONNECTED-LAKE INTERCHANGE VOLUMES
C
      FLSUM = ZERO
      DO 1325 L=1,NLAKES
      FLSUM = FLSUM + FLXINL(L)
 1325 CONTINUE
      IF(LWRT.GT.0.OR.ICBCFL.LE.0) GO TO 1350
      TV = TVOLM/1000000.
      IF(FLSUM.GE.TV) WRITE(IOUT,1327) FLSUM
 1350 CONTINUE
C
C   CHECK FOR LAKE CELLS GOING DRY
C
      LDR = 0
      DO 950 L=1,LKNODE
         IL=ILAKE(1,L)
         IR=ILAKE(2,L)
         IC=ILAKE(3,L)
         LAKE=ILAKE(4,L)
C  Convert ILAKE(5,L):  1 and 2 are type 1,  3 and 4 are type 2, 6 is type 0
         ITYPE = (ILAKE(5,L)+1)/2
         IF(ITYPE.EQ.3) ITYPE=0
         IF(ITYPE.NE.0) GO TO 950
         IF(IBOUND(IC,IR,IL).GT.0) BOTLK = BOTM(IC,IR,LBOTM(IL-1))
         IF(IBOUND(IC,IR,IL).EQ.0) BOTLK = BOTM(IC,IR,LBOTM(IL))
         IF (STGNEW(LAKE).LE.BOTLK) LDR = LDR + 1
         IF (STGNEW(LAKE).LE.BOTLK) LDRY(LDR) = L
  950 CONTINUE
      IF(LWRT.GT.0.OR.ICBCFL.LE.0) GO TO 951
      IF(LDR.EQ.0) WRITE(IOUT,875) LDR
  875 FORMAT(//1X,I5,2X,'LAKE CELLS ARE DRY.')
      IF(LDR.EQ.0) GO TO 951
      WRITE(IOUT,874) LDR
  874 FORMAT(/5X,'SECTIONS OF THE LAKE BOTTOM HAVE BECOME DRY.  THE DRY
     1SECTIONS'/5X,'LIE ABOVE THE FOLLOWING',I3,' AQUIFER CELLS (LAYER,R
     2OW,COLUMN):')
      LDR1 = 0
      DO 952 L=1,LDR
      LDR1 = LDR1 + 1
      L1 = LDRY(L)
      ILB(LDR1) = ILAKE(1,L1)
      IRB(LDR1) = ILAKE(2,L1)
      ICB(LDR1) = ILAKE(3,L1)
      IF(LDR1.LT.5) GO TO 952
      WRITE(IOUT,876) (ILB(I),IRB(I),ICB(I),I=1,5)
  876 FORMAT(5X,5('(',I3,',',I3,',',I3,')',2X))
      LDR1 = 0
  952 CONTINUE
      IF(LDR1.GT.0) WRITE(IOUT,876) (ILB(I),IRB(I),ICB(I),I=1,LDR1)
  951 CONTINUE
      IF(IUNITGWT.GT.0) GO TO 1086
      NDUM=1
      IF(IUNITGAGE.GT.0) CALL GAGE5LO(IGGLST,NUMGAGE,IUNITGWT,STGNEW,
     *  BUFF,NLAKES,TOTIM,NDUM,VOL,
     *                  PRECIP,EVAP,RNF,
     *                  GWIN,GWOUT,SURFIN,SURFOT,
     *                  WTHDRW,FLXINL,SUMCNN,
     *                  STGOLD2,VOLOLD,STAGES,VOLINIT,BUFF,BUFF)
 1086 IF(LWRT.GT.0.OR.ICBCFL.LE.0) GO TO 1061
C
C   WRITE BUDGET SUMMARIES.
C
            WRITE(IOUT,1025) KPER, KSTP, DELT, PERTIM, TOTIM
 1025 FORMAT(/1X,'PERIOD ',I5,5X,'TIME STEP ',I5,5X,'TIME STEP LENGTH ',
     1   1PE11.4/1X,'PERIOD TIME ',E11.4,5X,'TOTAL SIMULATION TIME ',
     2   E11.4)
            WRITE(IOUT,1040)
 1040 FORMAT(//17X,'HYDROLOGIC BUDGET SUMMARIES FOR SIMULATED LAKES'
     1    ,/,5X,'(ALL FLUID FLUXES ARE VOLUMES ADDED TO THE LAKE DURING'
     2    ,' PRESENT TIME STEP)'
     3    ,/,5X,'------------------------------------------------------'
     4    ,'-------------------')
            WRITE(IOUT,1045)
 1045 FORMAT(1X,'LAKE',3X,'STAGE',5X,'PRECIP',7X,'EVAP',7X,'RUNOFF')
 1061       DO 1100 NN=1,NLAKES
               PPTIN=PRECIP(NN)
               EOUT=EVAP(NN)
       IF(RNF(NN).GE.ZERO) RUNF = RNF(NN)
       IF(RNF(NN).LT.ZERO) RUNF =-RNF(NN)*PRCPLK(NN)*BGAREA(NN)
         RUNFD = RUNF*DELT
C
                  CUMPPT(NN)=CUMPPT(NN)+PPTIN
                  CUMEVP(NN)=CUMEVP(NN)+EOUT
                  CUMRNF(NN)=CUMRNF(NN)+RUNFD
                  IF(LWRT.GT.0.OR.ICBCFL.LE.0) GO TO 1100
                  WRITE(IOUT,1050) NN,STGNEW(NN),PPTIN,EOUT,RUNFD
 1100       CONTINUE
 1050 FORMAT(1X,I3,1X,F8.2,1P,3E12.4)
       IF(LWRT.LE.0.AND.ICBCFL.GT.0) WRITE(IOUT,1046)
 1046 FORMAT(/12X,'GROUND WATER',12X,'SURFACE WATER'/1X,'LAKE',4X,
     1 'INFLOW',5X,'OUTFLOW',6X,'INFLOW',5X,'OUTFLOW')
            DO 1101 NN=1,NLAKES
               QIN=GWIN(NN)*DELT
               QOUT=GWOUT(NN)*DELT
               QSIN=SURFIN(NN)
               QSOUT=SURFOT(NN)
C
                   CUMGWI(NN)=CUMGWI(NN)+QIN
                   CUMGWO(NN)=CUMGWO(NN)+QOUT
                   CUMSWI(NN)=CUMSWI(NN)+QSIN
                   CUMSWO(NN)=CUMSWO(NN)+QSOUT
 1101              IF(LWRT.LE.0.AND.ICBCFL.GT.0)
     1                 WRITE(IOUT,1051) NN,QIN,QOUT,QSIN,QSOUT
 1051 FORMAT(1X,I3,1P,4E12.4)
      IF(LWRT.LE.0.AND.ICBCFL.GT.0) WRITE(IOUT,1035)
 1035 FORMAT(/9X,'WATER',4X,'CONNECTED LAKE',3X,'UPDATED',5X,'TIME-STEP'
     1 ,10X,'STAGE CHANGE'/1X,'LAKE',5X,'USE',9X,'INFLUX',
     2  8X,'VOLUME',4X,'SURFACE AREA',3X,'TIME STEP',2X,'CUMULATIVE')
      DO 1105 NN=1,NLAKES
        WDRAW=WTHDRW(NN)*DELT
C
        CUMWDR(NN)=CUMWDR(NN)+WDRAW
        CUMFLX(NN)=CUMFLX(NN)+FLXINL(NN)
 1105 IF(LWRT.LE.0.AND.ICBCFL.GT.0)
     1       WRITE(IOUT,1055) NN,WDRAW,FLXINL(NN),VOL(NN),SURFA(NN),
     2           DELH(NN),TDELH(NN)
 1055 FORMAT(1X,I3,1P,3E13.4,1X,2E13.4,E12.4)
      IF(LWRT.GT.0.OR.ICBCFL.LE.0) GO TO 1041
      WRITE(IOUT,1301)
C
C   WRITE CUMULATIVE BUDGET SUMMARIES.
C
            WRITE(IOUT,2040)
 2040 FORMAT(//12X,'CUMULATIVE HYDROLOGIC BUDGET SUMMARIES FOR SIMULATED
     1 LAKES'
     2    ,/,5X,'(ALL FLUID FLUXES ARE SUMS OF VOLUMES ADDED TO THE'
     3    ,' LAKE SINCE INITIAL TIME)'
     4    ,/,5X,'------------------------------------------------------'
     5    ,'---------------------')
            WRITE(IOUT,2045)
 2045 FORMAT(1X,'LAKE',7X,'PRECIP',7X,'EVAP',7X,'RUNOFF')
            DO 2100 NN=1,NLAKES
 2100 WRITE(IOUT,2050) NN,CUMPPT(NN),CUMEVP(NN),CUMRNF(NN)
 2050 FORMAT(1X,I3,3X,1P,3E12.4)
            WRITE(IOUT,2046)
 2046 FORMAT(/12X,'GROUND WATER',12X,'SURFACE WATER'/1X,'LAKE',4X,
     1 'INFLOW',5X,'OUTFLOW',6X,'INFLOW',5X,'OUTFLOW')
            DO 2101 NN=1,NLAKES
 2101 WRITE(IOUT,2051) NN,CUMGWI(NN),CUMGWO(NN),CUMSWI(NN),CUMSWO(NN)
 2051 FORMAT(1X,I3,1P,4E12.4)
      WRITE(IOUT,2035)
 2035 FORMAT(/9X,'WATER',4X,'CONNECTED LAKE'/
     1       1X,'LAKE',5X,'USE',9X,'INFLUX')
      DO 2105 NN=1,NLAKES
 2105 WRITE(IOUT,2055) NN,CUMWDR(NN),CUMFLX(NN)
 2055 FORMAT(1X,I3,1P,2E13.4)
      WRITE(IOUT,1301)
      IF(KCNT.LE.0) GO TO 11041
      IF (KCNT.GT.1) THEN
        WRITE(IOUT,11055) KCNT
11055 FORMAT(/1X,I3,' CONNECTED LAKE SETS'/)
        DO 11056 IIC=1,KCNT
          JIC = MSUB1(IIC)
          WRITE(IOUT,11057) JIC, (MSUB(LIC,IIC),LIC=1,JIC)
11057 FORMAT(1X,I3,' LAKES:  ',25I3)
11056   CONTINUE
      ELSE
        WRITE(IOUT,21055) KCNT
21055 FORMAT(/1X,I3,' CONNECTED LAKE SET'/)
        JIC = MSUB1(IIC)
        WRITE(IOUT,11057) JIC, (MSUB(LIC,IIC),LIC=1,JIC)
      END IF
11041 CONTINUE
 1041 CONTINUE
C
C13-----IF C-B-C TERMS WILL BE SAVED THEN WRITE TO DISK.
      IF(IBD.EQ.1) CALL UBUDSV(KSTP,KPER,TEXT,ILKCB,BUFF,NCOL,NROW,
     1                         NLAY,IOUT)
C
C14A-----MOVE RATES INTO VBVL FOR PRINTING BY MODULE BAS OT.
1200  VBVL(3,MSUM)=RATIN
      VBVL(4,MSUM)=RATOUT
C
C14B-----MOVE PRODUCT OF RATE AND TIME STEP INTO VBVL ACCUMULATORS.
      VBVL(1,MSUM)=VBVL(1,MSUM)+RATIN*DELT
      VBVL(2,MSUM)=VBVL(2,MSUM)+RATOUT*DELT
C
C14C-----MOVE BUDGET TERM LABELS INTO VBVM FOR PRINTING BY BAS OT.
      VBNM(MSUM)=TEXT
C15------INCREASE BUDGET COUNTER.
      MSUM=MSUM+1
C
C12-----RETURN.
      RETURN
      END
C
      SUBROUTINE LSN3RP(NTRB,NDV,NLAKES,ITRB,IDIV,NSS,IDIVAR,IOTSG,
     1                  IOUT,NODES,RK)
C
C    *******************************************************************
C--  IF STREAMS EXIST, DEFINE CONNECTIONS BETWEEN LAKES AND STREAMS
C    *******************************************************************
C
C    -------------------------------------------------------------------
C        SPECIFICATIONS:
C
      DIMENSION ITRB(NLAKES,NSS),IDIV(NLAKES,NSS),IOTSG(NSS)
      DIMENSION RK(2,NLAKES),IDIVAR(2,NSS)
C    -------------------------------------------------------------------
C
C-- DOUBLE CHECK SIZE OF RK (STORED IN BUFF) vs. NLAKES
C
      IF ((NLAKES*2).GT.NODES) THEN
         WRITE (IOUT,*) '***NLAKES too large for BUFF in Subroutine LSN3
     1RP***  STOP EXECUTION'
        call stopfile ! emrl
         STOP
      END IF
C
C-- INITIALIZE ARRAYS
C
      DO 50 I=1,NSS
      DO 50 LK=1,NLAKES
      ITRB(LK,I) = 0
   50 IDIV(LK,I) = 0
      DO 55 LK=1,NLAKES
      RK(1,LK) = 0.
   55 RK(2,LK) = 0.
      NTRB = 0
      NDV = 0
C
C-- Build arrays to define lake tributary & diversion links ...
C        based on stream package input data
C
C---  Stream Inflow to Lakes
      DO 100 ISEG=1,NSS
      IF(IOTSG(ISEG).LT.0) THEN
        LAKE = -IOTSG(ISEG)
        RK(1,LAKE) = RK(1,LAKE) + 1.
        K1 = RK(1,LAKE)
        ITRB(LAKE,K1) = ISEG
        IF(RK(1,LAKE).GT.NTRB) NTRB = RK(1,LAKE)
      ENDIF
C
C---  Stream Outflow from Lakes
      IF(IDIVAR(1,ISEG).LT.0) THEN
        LAKE = -IDIVAR(1,ISEG)
        RK(2,LAKE) = RK(2,LAKE) + 1.
        K1 = RK(2,LAKE)
        IDIV(LAKE,K1) = ISEG
        IF(RK(2,LAKE).GT.NDV) NDV = RK(2,LAKE)
      ENDIF
  100 CONTINUE
C
C--  PRINT LAKE INFLOW STREAM SEGMENTS.
      WRITE(IOUT,10)
10    FORMAT(6X,'LAKE ',4X,'INFLOWING STREAM SEGMENT')
      DO 520 IK=1,NLAKES
      DO 519 JK=1,NSS
      IF(ITRB(IK,JK).LE.0) GO TO 521
  519 CONTINUE
  521 JK1 = JK - 1
      IF(JK1.GT.0) WRITE(IOUT,15) IK,(ITRB(IK,JK),JK=1,JK1)
15    FORMAT(5X,I5,14X,100I5)
  520 CONTINUE
      WRITE(IOUT,103) NTRB
103    FORMAT(/1X,'MAXIMUM NUMBER OF STREAMS INFLOWING TO A',
     1    ' LAKE IS',I5/)
C
C--  PRINT LAKE STREAM OUTFLOW SEGMENT (FROM A LAKE) NUMBERS.
C
      WRITE(IOUT,13)
13    FORMAT(6X,'LAKE ',4X,'OUTFLOWING STREAM',' SEGMENT')
      DO 600 IK=1,NLAKES
      DO 523 JK=1,NSS
      IF(IDIV(IK,JK).LE.0) GO TO 527
  523 CONTINUE
  527 JK1 = JK - 1
      IF(JK1.GT.0) WRITE(IOUT,15) IK,(IDIV(IK,JK),JK=1,JK1)
  600 CONTINUE
      WRITE(IOUT,133) NDV
133   FORMAT(/1X,'MAXIMUM NUMBER OF STREAMS OUTFLOWING',
     1    ' FROM A LAKE IS',I5/)
C
C-- RETURN
      RETURN
      END
      SUBROUTINE GWF1LAK3BCF6RP(IOUT,BEDLAK,LKARR1,ILAKE,CNDFCT,
     1         DELR,DELC,HY,TRPY,LAYHDT,MXLKND,NCOL,NROW,NLAY,
     2         LKNODE,IWDFLG,CVWD)
C
C     ******************************************************************
C     COMPUTE VERTICAL CONDUCTANCES AND HORIZONTAL CONDUCTANCES PER UNIT
C     THICKNESS FOR LAKES WHEN BCF PACKAGE IS USED
C     ******************************************************************
C
C     ------------------------------------------------------------------
C     SPECIFICATIONS:
      DIMENSION CNDFCT(MXLKND),BEDLAK(MXLKND),ILAKE(5,MXLKND)
      DIMENSION LKARR1(NCOL,NROW,NLAY),HY(NCOL,NROW,NLAY),
     1          CVWD(NCOL,NROW,NLAY)
      DIMENSION DELR(NCOL),DELC(NROW),TRPY(NLAY),LAYHDT(NLAY)
C     ------------------------------------------------------------------
C
      WRITE(IOUT,108)
  108 FORMAT(//9X,'C',15X,'INTERFACE CONDUCTANCES BETWEEN LAKE AND ',
     1  'AQUIFER CELLS'/
     2  3X,'L',5X,'O',10X,'(IF TYPE = 0, CONDUCTANCE (L^2/T) IS ',
     3  'BETWEEN AQUIFER CELL AND OVERLYING LAKE CELL.)',/
     4  3X,'A',5X,'L',2X,'L',2X,'T',
     5  4X,'(IF TYPE = 1 TO 4, CONDUCTANCES ARE PER UNIT SATURATED ',
     6  'THICKNESS (L/T).)'/
     7  3X,'Y',2X,'R',2X,'U',2X,'A',2X,'Y'/
     8  3X,'E',2X,'O',2X,'M',2X,'K',2X,'P',
     9  24X,'LAKEBED',6X,'C O N D U C T A N C E S'/3X,'R',2X,'W',2X,
     1  'N',2X,'E',
     2  2X,'E',5X,'DELTA Y',3X,'DELTA X',2X,'LEAKANCE',3X,'LAKEBED',3X,
     3  'AQUIFER',2X,'COMBINED'/1X,79('_'))
C
      DO 350 II=1,LKNODE
      K = ILAKE(1,II)
      J = ILAKE(2,II)
      I = ILAKE(3,II)
      CNDFCT(II) = 0.0
C  Convert ILAKE(5,II):  1 and 2 are type 1,  3 and 4 are type 2, 6 is type 0
      NTYP = (ILAKE(5,II)+1)/2
      IF(NTYP.EQ.3) NTYP=0
      NTYP = NTYP + 1
      IF(NTYP.EQ.1) THEN
C
C  Vertical Conductance
C    for vertical interface, "K" is layer below bottom of lake
        CNDFC1=0.0
        IF(K.EQ.NLAY.AND.LKARR1(I,J,K).GT.0) GO TO 315
        IF(BEDLAK(II).LE.0.0) GO TO 315
        CNDFC1 = BEDLAK(II)*DELR(I)*DELC(J)
        IF (IWDFLG.EQ.0) THEN
          CNDFCT(II) = CNDFC1
        ELSE
          IF(CVWD(I,J,K-1).LE.0.0.OR.CNDFC1.LE.0.0) GO TO 350
          CNDFCT(II) = 1.0/(0.5/CVWD(I,J,K-1)+1.0/CNDFC1)
        END IF
  315   IF (IWDFLG.EQ.0) THEN
          WRITE(IOUT,7324) (ILAKE(I1,II),I1=1,5),DELC(J),DELR(I),
     1        BEDLAK(II),CNDFC1,CNDFCT(II)
 7324     FORMAT(1X,5I3,2X,1P,4E10.2,10X,E10.2)
        ELSE
          CVWD2= 2.0*CVWD(I,J,K-1)
          WRITE(IOUT,7325) (ILAKE(I1,II),I1=1,5),DELC(J),DELR(I),
     1        BEDLAK(II),CNDFC1,CVWD2,CNDFCT(II)
 7325     FORMAT(1X,5I3,2X,1P,6E10.2)
        END IF
      ELSE
C
C  Horizontal conductance
        IF(LAYHDT(K).EQ.0) THEN
          WRITE(IOUT,346) I,J,K
  346     FORMAT(1X,
     1    'ERROR -- NODE ADJACENT TO LAKE IS IN CONFINED LAYER',5X,3I5)
        call stopfile ! emrl
          STOP
        END IF
C
        TT = HY(I,J,K)
        IF(NTYP.EQ.2) CNDFC2 = 2.0*TT*DELC(J)/DELR(I)
        IF(NTYP.EQ.3) CNDFC2 = 2.0*TRPY(K)*TT*DELR(I)/DELC(J)
        IF(NTYP.EQ.2) CNDFC1 = BEDLAK(II)*DELC(J)
        IF(NTYP.EQ.3) CNDFC1 = BEDLAK(II)*DELR(I)
        CNDFCT(II) = 1.0/(1.0/CNDFC2+1.0/CNDFC1)
        WRITE(IOUT,7325) (ILAKE(I1,II),I1=1,5),DELC(J),DELR(I),
     1    BEDLAK(II),CNDFC1,CNDFC2,CNDFCT(II)
      END IF
  350 CONTINUE
C
      RETURN
      END
C
      SUBROUTINE GWF1LAK3LPF1RP(IOUT,BEDLAK,LKARR1,ILAKE,CNDFCT,
     1         DELR,DELC,HK,HANI,LAYHDT,MXLKND,NCOL,NROW,NLAY,
     2         LKNODE,VKA,VKCB,BOTM,NBOTM)
C
C     ******************************************************************
C     COMPUTE VERTICAL CONDUCTANCES AND HORIZONTAL CONDUCTANCES PER UNIT
C     THICKNESS FOR LAKES WHEN LPF PACKAGE IS USED
C     ******************************************************************
C
C     ------------------------------------------------------------------
C     SPECIFICATIONS:
      DIMENSION CNDFCT(MXLKND),BEDLAK(MXLKND),ILAKE(5,MXLKND)
      DIMENSION LKARR1(NCOL,NROW,NLAY),HK(NCOL,NROW,NLAY),
     1          VKA(NCOL,NROW,NLAY),BOTM(NCOL,NROW,0:NBOTM),
     2          VKCB(NCOL,NROW,NLAY)
      DIMENSION DELR(NCOL),DELC(NROW),HANI(NCOL,NROW,NLAY),
     1          LAYHDT(NLAY)
C
      COMMON /DISCOM/LBOTM(200),LAYCBD(200)
      COMMON /LPFCOM/LAYTYP(200),LAYAVG(200),CHANI(200),LAYVKA(200),
     1               LAYWET(200)
C     ------------------------------------------------------------------
C
      WRITE(IOUT,108)
  108 FORMAT(//9X,'C',15X,'INTERFACE CONDUCTANCES BETWEEN LAKE AND ',
     1  'AQUIFER CELLS'/
     2  3X,'L',5X,'O',10X,'(IF TYPE = 0, CONDUCTANCE (L^2/T) IS ',
     3  'BETWEEN AQUIFER CELL AND OVERLYING LAKE CELL.)',/
     4  3X,'A',5X,'L',2X,'L',2X,'T',
     5  4X,'(IF TYPE = 1 TO 4, CONDUCTANCES ARE PER UNIT SATURATED ',
     6  'THICKNESS (L/T).)'/
     7  3X,'Y',2X,'R',2X,'U',2X,'A',2X,'Y'/
     8  3X,'E',2X,'O',2X,'M',2X,'K',2X,'P',
     9  24X,'LAKEBED',6X,'C O N D U C T A N C E S'/3X,'R',2X,'W',2X,
     1  'N',2X,'E',
     2  2X,'E',5X,'DELTA Y',3X,'DELTA X',2X,'LEAKANCE',3X,'LAKEBED',3X,
     3  'AQUIFER',2X,'COMBINED'/1X,79('_'))
C
      DO 350 II=1,LKNODE
      K = ILAKE(1,II)
      J = ILAKE(2,II)
      I = ILAKE(3,II)
      CAQ = 0.0
      CNDFCT(II) = 0.0
C  Convert ILAKE(5,II):  1 and 2 are type 1,  3 and 4 are type 2, 6 is type 0
      NTYP = (ILAKE(5,II)+1)/2
      IF(NTYP.EQ.3) NTYP=0
      NTYP=NTYP + 1
      IF(NTYP.EQ.1) THEN
C
C  Vertical Conductance
C    for vertical interface, "K" is layer below bottom of lake
        CNDFC1=0.0
        IF(K.EQ.NLAY.AND.LKARR1(I,J,K).GT.0) GO TO 315
        IF(BEDLAK(II).LE.0.0) GO TO 315
        CNDFC1 = BEDLAK(II)*DELR(I)*DELC(J)
        IF(LAYVKA(K).EQ.0) THEN
           VK=VKA(I,J,K)
        ELSE
           VK=HK(I,J,K)/VKA(I,J,K)
        END IF
c   skip if zero vk
        IF(VK.LE.0.0) GO TO 350
        BBOT=BOTM(I,J,LBOTM(K))
        TTOP=BOTM(I,J,LBOTM(K)-1)
        CAQ=VK*DELR(I)*DELC(J)/((TTOP-BBOT)*0.5)
        IF(LAYCBD(K-1).GT.0) THEN
c   skip if zero vkcb
          IF(VKCB(I,J,LAYCBD(K)).LE.0.0) GO TO 350
          BBOT=BOTM(I,J,LBOTM(K)-1)
          TTOP=BOTM(I,J,LBOTM(K-1))
          CCB=VKCB(I,J,LAYCBD(K-1))*DELR(I)*DELC(J)/(TTOP-BBOT)
          !include VKCB
          CAQ = 1.0/(1.0/CAQ + 1.0/CCB)
        END IF
        CNDFCT(II) = 1.0/(1.0/CAQ+1.0/CNDFC1)
  315   WRITE(IOUT,7325) (ILAKE(I1,II),I1=1,5),DELC(J),DELR(I),
     1             BEDLAK(II),CNDFC1,CAQ,CNDFCT(II)
      ELSE
C
C  Horizontal conductance
        IF(LAYHDT(K).EQ.0) THEN
          WRITE(IOUT,346) I,J,K
  346     FORMAT(1X,
     1    'ERROR -- NODE ADJACENT TO LAKE IS IN CONFINED LAYER',5X,3I5)
        call stopfile ! emrl
          STOP
        END IF
C
        TT = HK(I,J,K)
C X-DIRECTION
        IF(NTYP.EQ.2) CNDFC2 = 2.0*TT*DELC(J)/DELR(I)
C Y-DIRECTION
        IF(NTYP.EQ.3) THEN
          IF(CHANI(K).LE.0) THEN
            KH=-CHANI(K)
            CNDFC2 = 2.0*HANI(I,J,KH)*TT*DELR(I)/DELC(J)
          ELSE
            CNDFC2 = 2.0*CHANI(K)*TT*DELR(I)/DELC(J)
          END IF
        END IF
        IF(NTYP.EQ.2) CNDFC1 = BEDLAK(II)*DELC(J)
        IF(NTYP.EQ.3) CNDFC1 = BEDLAK(II)*DELR(I)
        CNDFCT(II) = 1.0/(1.0/CNDFC2+1.0/CNDFC1)
        WRITE(IOUT,7325) (ILAKE(I1,II),I1=1,5),DELC(J),DELR(I),
     1    BEDLAK(II),CNDFC1,CNDFC2,CNDFCT(II)
 7325   FORMAT(1X,5I3,2X,1P,6E10.2)
      END IF
  350 CONTINUE
C
      RETURN
      END
C
      SUBROUTINE GWF1LAK3HUF1RP(IOUT,BEDLAK,LKARR1,ILAKE,CNDFCT,
     1         DELR,DELC,HK,HKCC,LAYHDT,MXLKND,NCOL,NROW,NLAY,
     2         LKNODE,VKA,BOTM,NBOTM)
C
C     ******************************************************************
C     COMPUTE VERTICAL CONDUCTANCES AND HORIZONTAL CONDUCTANCES PER UNIT
C     THICKNESS FOR LAKES WHEN HUF PACKAGE IS USED
C     ******************************************************************
C
C     ------------------------------------------------------------------
C     SPECIFICATIONS:
      DIMENSION CNDFCT(MXLKND),BEDLAK(MXLKND),ILAKE(5,MXLKND)
      DIMENSION LKARR1(NCOL,NROW,NLAY),HK(NCOL,NROW,NLAY),
     1          VKA(NCOL,NROW,NLAY),BOTM(NCOL,NROW,0:NBOTM)
      DIMENSION DELR(NCOL),DELC(NROW),HKCC(NCOL,NROW,NLAY),
     1          LAYHDT(NLAY)
C
      COMMON /DISCOM/LBOTM(200),LAYCBD(200)
      COMMON /LPFCOM/LAYTYP(200),LAYAVG(200),CHANI(200),LAYVKA(200),
     1               LAYWET(200)
C     ------------------------------------------------------------------
C
      WRITE(IOUT,108)
  108 FORMAT(//9X,'C',15X,'INTERFACE CONDUCTANCES BETWEEN LAKE AND ',
     1  'AQUIFER CELLS'/
     2  3X,'L',5X,'O',10X,'(IF TYPE = 0, CONDUCTANCE (L^2/T) IS ',
     3  'BETWEEN AQUIFER CELL AND OVERLYING LAKE CELL.)',/
     4  3X,'A',5X,'L',2X,'L',2X,'T',
     5  4X,'(IF TYPE = 1 TO 4, CONDUCTANCES ARE PER UNIT SATURATED ',
     6  'THICKNESS (L/T).)'/
     7  3X,'Y',2X,'R',2X,'U',2X,'A',2X,'Y'/
     8  3X,'E',2X,'O',2X,'M',2X,'K',2X,'P',
     9  24X,'LAKEBED',6X,'C O N D U C T A N C E S'/3X,'R',2X,'W',2X,
     1  'N',2X,'E',
     2  2X,'E',5X,'DELTA Y',3X,'DELTA X',2X,'LEAKANCE',3X,'LAKEBED',3X,
     3  'AQUIFER',2X,'COMBINED'/1X,79('_'))
C
      DO 350 II=1,LKNODE
      K = ILAKE(1,II)
      J = ILAKE(2,II)
      I = ILAKE(3,II)
      CAQ = 0.0
      CNDFCT(II) = 0.0
C  Convert ILAKE(5,II):  1 and 2 are type 1,  3 and 4 are type 2, 6 is type 0
      NTYP = (ILAKE(5,II)+1)/2
      IF(NTYP.EQ.3) NTYP=0
      NTYP=NTYP + 1
      IF(NTYP.EQ.1) THEN
C
C  Vertical Conductance
C    for vertical interface, "K" is layer below bottom of lake
        CNDFC1=0.0
        IF(K.EQ.NLAY.AND.LKARR1(I,J,K).GT.0) GO TO 315
        IF(BEDLAK(II).LE.0.0) GO TO 315
        CNDFC1 = BEDLAK(II)*DELR(I)*DELC(J)
        VK=VKA(I,J,K)
c   skip if zero vk
        IF(VK.LE.0.0) GO TO 350
        BBOT=BOTM(I,J,LBOTM(K))
        TTOP=BOTM(I,J,LBOTM(K)-1)
        CAQ=VK*DELR(I)*DELC(J)/((TTOP-BBOT)*0.5)
        CNDFCT(II) = 1.0/(1.0/CAQ+1.0/CNDFC1)
  315   WRITE(IOUT,7325) (ILAKE(I1,II),I1=1,5),DELC(J),DELR(I),
     1             BEDLAK(II),CNDFC1,CAQ,CNDFCT(II)
      ELSE
C
C  Horizontal conductance
        IF(LAYHDT(K).EQ.0) THEN
          WRITE(IOUT,346) I,J,K
  346     FORMAT(1X,
     1    'ERROR -- NODE ADJACENT TO LAKE IS IN CONFINED LAYER',5X,3I5)
        call stopfile ! emrl
          STOP
        END IF
C
        TT = HK(I,J,K)
        TY = HKCC(I,J,K)
C X-DIRECTION
        IF(NTYP.EQ.2) CNDFC2 = 2.0*TT*DELC(J)/DELR(I)
C Y-DIRECTION
        IF(NTYP.EQ.3) CNDFC2 = 2.0*TY*DELC(J)/DELR(I)
        IF(NTYP.EQ.2) CNDFC1 = BEDLAK(II)*DELC(J)
        IF(NTYP.EQ.3) CNDFC1 = BEDLAK(II)*DELR(I)
        CNDFCT(II) = 1.0/(1.0/CNDFC2+1.0/CNDFC1)
        WRITE(IOUT,7325) (ILAKE(I1,II),I1=1,5),DELC(J),DELR(I),
     1    BEDLAK(II),CNDFC1,CNDFC2,CNDFCT(II)
 7325   FORMAT(1X,5I3,2X,1P,6E10.2)
      END IF
  350 CONTINUE
C
      RETURN
      END
